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Abstract 

Relativistic Heavy Ion Collisions: 
Viscous Hydrodynamic Simulations and Final State Interactions 

Matthew W. Luzum 

Chair of the Supervisory Committee: 
Professor Gerald A. Miller 
Department of Physics 

In this dissertation I introduce relativistic heavy ion collisions and describe theoretical 
approaches to understanding them — in particular, viscous hydrodynamic simulations and 
investigations of final state interactions. 

The successful ideal hydrodynamic models of the collisions at the Relativistic Heavy 
Ion Collider (RHIC) were extended by performing viscous hydrodynamic simulations. This 
was done by making use of the recently derived full conformally invariant second order 
relativistic viscous hydrodynamic equations. Results for multiplicity, radial flow and elliptic 
flow in ^sjsfN = 200 GeV Au+Au RHIC collisions are presented and the range of the ratio 
of shear viscosity over entropy density ^ for which our hydrodynamic model is consistent 
with experimental data is quoted. 

In addition, simulations were performed of the planned y^sjvjv = 5.5 TeV Pb+Pb and 
^/s = 14 TeV p+p collisions at the Large Hadron Collider (LHC). The elliptic flow coefficient 
V2 is predicted to be 10% larger for the Pb+Pb collisions compared to top energy RHIC 
collisions, and is predicted to be consistent with zero for proton collisions unless ^ < 0.08. 

Finally, final state interactions were investigated within the distorted wave emission 
function (DWEF) model. Work is presented on an improved understanding of the DWEF 
model, and the potential effect of final state interactions in the form of a pion optical 
potential on the elliptic flow coefficient V2 was calculated to be at the ~ 20% level. 
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Chapter 1 

PROLOGUE: INTRODUCTION, MOTIVATION AND PHILOSOPHY 

The structure of this dissertation is as fohows. This chapter consists of a nontechnical 
and general discussion of the motivation behind relativistic heavy ion collisions. Chapter 
[2] briefly describes the experiments and introduces a few measured quantities that will 
be important for the theoretical work that is presented in the remaining chapters. The 
framework of hydrodynamic theory in general is introduced in chapter 3} Following these 



introductory chapters is the main body, which presents original work (collaboratively) done 



by the author. In chapter 4, hydrodynamic models of collisions at the Relativistic Heavy 
Ion Coh 



Refs 



ider are introduced and results of these simulations are presented (corresponding to 

nh , , 

. [l|, y] ) . Simulations of collisions at the Large Hadron Collider are given in [chapter "5 



(corresponding to Ref. [3]). Chapters [6] and [7] describe the DWEF model investigations of 
final state interactions (roughly corresponding to Refs. P, 0]). AppendixlAl offers additional 
details of the DWEF calculation of f 2 while appendix [B] contains a list of the conventions 
and notation used throughout this manuscript. 

A digital version of this manuscript, including high quality color figures, will be available 



online at http : / / arxiv . org/a/luzuin_in_l At the time of this writing, source code and 



results from the viscous hydrodynamic simulations presented in chapters [H and [5] can be 



found on Paul Romatschke's webpage: http : / /hep . itp . tuwien . ac . at/^aulrom/ 
1.1 Introduction 

The goal of the work described in this dissertation is to better understand how the world 
works on the most fundamental levels. By studying very small and/or simple physical 
systems, we can extract information about the fundamental laws that govern the world we 
live in — including (presumably) the behavior of the much more complex systems that we 
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typically encounter in day-to-day life. Wc do this because of simple human curiosity and 
a natural desire for knowledge, but also because this knowledge tends to be very useful. 
When we have a detailed understanding of how the world works, we can often manipulate 
it for our benefit. 

Of course this particular line of scientific inquiry is only one of many that are both useful 
and necessary. More complex systems must be studied on their own as well. For example, 
it is neither interesting nor practically possible to calculate the fundamental interactions of 
every molecule in a bridge when trying to determine whether it will support a load without 
collapsing (let alone all the atoms and electrons or quarks and gluons). Even in the work 
described herein, hydrodynamic equations will be used extensively. These equations describe 
a sort of coarse-grained behavior on a scale that is large compared to the fundamental 
microscopic dynamics to which hydrodynamic behavior is largely insensitive. Indeed, many 
would argue that the study of larger scale and perhaps less fundamental behavior — e.g., 
chemistry, biology, materials science, medicine, etc. — is more important. Nevertheless, I 
would argue that it is still very important to study these fundamental laws of nature — even 
in such exotic regimes as extremely high temperature nuclear matter, and even if it doesn't 
seem to have any obvious practical applications. A hundred years ago, there was no reason 
to think that understanding the weird quantum mechanical behavior of tiny particles would 
be of any practical use. On the contrary, almost none of the current technology that we all 
take for granted — and that enable much of the progress in other sciences — would be possible 
without the insights gained from these seemingly esoteric studies. 

1.2 Strong Interactions, QCD Phase Diagram and the Quark-Gluon Plasma 

1.2.1 What do we know about the world? 

The world as we know it is made up of matter and the forces that interact with and hold 
this matter together. These interactions are usually classified into four known fundamental 
forces: gravity, electricity and magnetism (electromagnetism) , the weak nuclear force, and 
the strong nuclear force. These fundamental forces are listed here in order of increasing 
strength, and therefore also generically of increasing importance as one considers behavior 
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at smaller and smaller length scales. 

For example, gravity is important for describing the movement of large collections of 
matter that have essentially no net electric charge, such as planets moving through the 
solar system. If we want to study how atoms form into molecules and solids, on the other 
hand, gravity has essentially no effect because electromagnetic interactions are so much 
stronger and are much more important to the movement of electrically charged matter. 
Going further down in scale, the structure of nuclei inside atoms is dominated by the strong 
and weak nuclear forces. 

The goal here will be to study particular aspects of the strong nuclear force, often referred 
to by physicists simply as strong interactions. Correspondingly, it will be necessary to look 
at very very small length scales. This will — due to Heisenberg's uncertainty principle — 
involve studying the behavior of matter at very large energies that are obtained, perhaps 
unsurprisingly, by smashing things together in accelerators. 

1.2.2 The strong force and quantum chromodynamics 

The Standard Model of elementary particle physics is composed of well-tested quantum field 
theories that describe all of the fundamental forces except gravity. The strong interactions 
in particular are well described by a theory called Quantum Chromodynamics (QCD). QCD 
describes the interaction of fundamental fields called quarks and gluons. These quarks and 
gluons form the protons and neutrons that, along with electrons, form almost all of the 
matter we see around us. 

In any particular interaction between fundamental particles in a quantum field theory 
such as QCD, there is a characteristic energy scale, and the strength of the interaction 
(quantified by a "coupling" g) depends on this scale. QCD has an unusual property called 
"asymptotic freedom", which means that the strength of the interaction decreases as this 
energy scale increases, and vice versa. An intrinsic energy scale for the strong interactions 
is Aqcd, where the coupling becomes order one. At energies much larger than this, the 
coupling is small and one can usually use the familiar methods of perturbation theory to 
calculate various quantities in QCD. Most of the precision tests of QCD have been done in 
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this regime and it is relatively well understood. When there are energy scales in a particular 
problem that are near or below Aqcd (even if the energies are very large compared to atomic 
energy scales), things typically become much more difficult, and this will be a hindrance in 
the study of heavy ion collisions. 

1.2.3 Confinement, temperature, and phase transitions 

Related to asymptotic freedom (but on the other end of the energy spectrum) is the concept 
of color confinement, another property of the strong interactions. Particles that participate 
in strong interactions have what's called a "color" charge, analogous to electric charge for 
electromagnetic behavior (and completely unrelated to the color of visible light). In loose 
terms, confinement means that it is impossible to isolate a color-charged particle such as 
a quark. They are only found tightly bound together with other colored objects in overall 
color-neutral states. 

As an example, think of a color-neutral pair of a quark and an antiquark. Asymptotic 
freedom implies that if they arc very close together they interact very weakly. If one was to 
try to pull them apart, however, the attraction would become stronger and stronger such 
that it would in principle take an infinite amount of energy to completely separate them. In 
reality there would eventually be so much energy between them that more quark-antiquark 
pairs would be created out of the vacuum, and you would just be left with multiple color- 
neutral states instead of the one you started with. 

One can imagine, however, collecting together some strongly interacting matter and rais- 
ing the temperature. Asymptotic freedom implies that there exists some (extremely large) 
temperature at which the strong interactions would become so weak that individual quarks 
and gluons would no longer be confined. Thus, there is expected to exist a deconfinement 
phase transition, where the color-neutral hadrons would "melt" into a deconfined state of 
matter called the quark-gluon plasma (QGP). 

How might one go about studying such an exotic regime of physics? Such huge temper- 
atures are hard to come by in the current universe — the temperature required (a few trillion 
degrees) is perhaps 100000 times larger than that at the center of the sun. The idea, then. 
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is to actually create conditions in the laboratory sufficient to create a QGP. By accelerating 
heavy ions (nuclei of large atoms) to extremely high energies and letting them collide, it 
was hoped that — if only for a tiny fraction of a second — the "fireball" created would be hot 
and dense enough to create this deconfined phase. By carefully studying what comes out 
of such a collision, one could then in principle discern many properties of such a state of 
matter and of the strong interactions in general. 
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Chapter 2 

RELATIVISTIC HEAVY ION COLLISIONS 
2.1 AGS, SPS, RHIC and LHC 

Once it was realized that QCD imphes a deconfined state of matter, the idea of colhding 
heavy ions to study it was quickly adopted. There already were relativistic heavy ion 
collisions being studied at the BEVALAC at Berkeley, but these were not at a sufficiently 
high energy to potentially create a QGP. The first ultrarelativistic heavy ion collisions that 
were done to investigate the quark-gluon plasma were performed at existing particle colliders 
that were modified to accept heavy ion beams, most notably the Alternating Gradient 
Synchrotron (AGS) at Brookhaven and the Super Proton Synchrotron (SPS) at CERN [6]. 
These were both fixed-target experiments that, among other things, collided Au+Au at up 
to 11 GeV per nucleon beam energy (AGS) and Pb+Pb at up to 160 A GeV (SPS). These 
experiments revealed tantalizing evidence of a hot and dense state of matter that had not 
previously been seen [3]. 

The first facility specifically designed for colliding ultrarelativistic heavy ions was the 
Relativistic Heavy Ion Collider (RHIC) at Brookhaven National Laboratory. There, colli- 
sions were performed by colliding beams with a center of mass energy of up to 200 GeV 
per nucleon pair, leading to much more energy potentially being deposited in the collision 
region as well as the possibility of the resulting fireball remaining in the deconfined state 
for a longer period of time, compared to the earlier lower energy collisions. All the experi- 
mental data analyzed in this dissertation are from runs at RHIC, and for simplicity the rest 
of chapter will focus on what has been done there. 

The Large Hadron Collider (LHC) at CERN is the next generation collider. Although 
most people know of it as a proton-proton collider, it was also designed to run heavy ion 
collisions part-time and the first runs should begin before long. Ultimately it is planned to 
collide Pb+Pb beams at up to 5.5 TeV per nucleon pair. 
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Figure 2.1: Thousands of particles are produced in a typical collision at RHIC. Picture is 
of particle tracks in a y^=200 A GeV Au+Au collision event from STAR collaboration [8]. 



2.2 Some Relevant Observables 



All the information that can be obtained from a given collision event comes from studying 
the thousands of produced particles that emerge from the collision region. Any information 
about the evolution of the fireball system and its medium properties must be inferred by 
looking at these final products well after they have stopped interacting with each other, as 



they stream into one of the detectors surrounding the collision region (see Figure 2.1 ). Much 
of what a theorist would ideally like to measure, therefore, may be inaccessible to direct 
measurement. However, a surprising amount about the collision can still be learned this 
way (see Refs. jg . [lol. Illl. [l^ for an overview from each of the main detector collaborations 



at RHIC of the first four years of results). The following describes the particular measured 
quantities that will be important for the work comprising this dissertation. 
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Figure 2.2: Cartoon of a heavy ion collision. The beam direction is the z direction (out of 
the page) while the x direction is in principle defined by the impact parameter h [left side 
of figure]. Anisotropic pressure gradients in the material deposited in a non-central (5 ^ 0) 
collision can lead to elliptic flow [right]. 



2.2.1 Single particle spectra and elliptic flow 

First, note some standard definitions: The beam direction defines the z-axis, and the x 



direction is defined such that the x-z plane is the collision plane as seen in Figure 2.2 It 
is often useful to use polar coordinates in the transverse (x-y) plane, where the azimuthal 
angle 4> is measured with respect to the x-axis. One can also define an impact parameter b 
connecting the lines-of-flight of the centers of mass of the colliding nuclei. 

The detected particles are characterized by their momentum p after they exit the collision 
region. Instead of the longitudinal component of momentum pz or the polar angle with 
respect to the beam 9, what is more commonly reported is the particle rapidity Y = 
arctanh{pz/E) or pseudorapidity r] = — ln[tan(|)]. (Note that rapidity and pseudorapidity 



are equivalent for a relativistically moving particle (m 0) since Y 



In 



and 



V 



iln 



Thus they are often used interchangeably, although one should always 



keep in mind the limits to the validity of this equivalence. 
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The probability of detecting a particle in a given event with a given rapidity and trans- 
verse momentum is given by the distribution ^^/'^ip^ = -^"^ such that 



where N is the total number of particles in the event and pT = yPx + Py is the transverse 
momentum. This distribution can refer to a particular species of identified particle (e.g., 
protons) or a composite measurement such as all charged hadrons combined. From this, 
one can calculate the mean transverse momentum 

It is useful to break up the distribution into its Fourier components with respect to the 
azimuthal angle of the outgoing particle's momentum (j)p = tan~"'^(^): 



dN 



dY(PpT 



1 + 2u„, cos(n 0p) . (2-3) 

n=l 

The sine terms vanish due to the reflection symmetry about the collision [x-z) plane, and 
for collisions of identical nuclei the reflection symmetry about the y-z plane causes the 
odd cosine moments to vanish. vq{pt) is referred to as radial flow and the next lowest 
non- vanishing coefficient V2{pt) is called the elliptic flow coefficient. Explicitly we have 

...(c„s(2«) = l/^c„.(2«-^. (2.5) 

The momentum integrated elliptic flow coefficient is denoted 

.r^ ^ff (2.6) 

and the minimum bias is defined by averaging over all collisions in a given run — i.e., 
integrating over impact parameter [3] 

b _ Jdbbv2{b) vo{b) 
" fdbbvoib) ■ ^^-^^ 
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One should note that the colliding nuclei consist of a collection of well-localized nucleons 



and so are not as smoothly distributed as Figure 2T2| might suggest. Therefore there can be 



some ambiguity in defining, e.g., the collision plane. (Also note that the odd moments of 
the particle distribution are not exactly zero). The theoretical calculations presented here, 
however, will model the collisions with smooth initial conditions and these definitions are 
then completely unambiguous. The difficulty then comes in knowing which experimental 



results to compare the theoretical results to. (See, e.g., chapter 4, comparing minimum bias 



V2 results to different experimental extractions that attempt to remove "non-ffow" effects). 

2.2.2 Two-particle correlations — Hanbury Brown/Twiss interferonietry 

In the 1950's, Robert Hanbury Brown and Richard Q. Twiss began using a method of 
intensity interferometry to measure the sizes of various objects — most notably measuring 
the angular size of the star Sirius in 1956 Q, [l^. It turns out that by correlating the 
intensity of light emitted incoherently from a source — even without measuring any infor- 
mation about phase (a mandatory ingredient of typical amplitude interferometry) — one can 
directly measure information about its size, as well the time dependence of a source that 
varies in time. This effect can be thought of as being caused by the symmetrization of the 
wavefunction for identical bosons, such as the photons emitted from a star, which causes 
enhancement of coincident measurement of pairs of these bosons with small momentum 
difference. (Similarly, an "anti-bunching" effect is present for identical fermions because of 
the antisymmetric nature of their wavefunction) . 

Despite early skepticism, this technique was quickly employed to measure space-time 



properties of collision systems in the laboratory 



identical particles emitted 
study heavy ion collisions 



:"rom such collisions 



3y analyzing two-particle correlations of 



16l |. and has since been used extensively to 



Such analyses of two-particle correlations are often generi- 
cally referred to as Hanbury Brown/Twiss (HBT) interferometry, or simply femtoscopy. 

Explicitly, the quantity constructed to analyze RHIC events is the ratio of the two- 
particle inclusive and single-particle inclusive spectra: 

_ dN/id^p^d^p2) 

^^P'^^'^= {dN/d^p,){dN/d^p,y ^'-'^ 
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The numerator measures the probabiUty that two particles of momentum pi and p2 are 
detected in the same event, while the denominator is the product of the familiar spectra 
from [section 2.2TT1 

Define the average momentum K = {pi +p2)/2 and the momentum difference q = 
(pi —P2)- Then one can define the directions L (longitudinal), O (out), and S (side) as the 
directions parallel to the beam, parallel to Kt = \/ + perpendicular to both 

the beam and Kt, respectively. 

The correlation function can then be parameterized as a Gaussian with parameters that 
are fit to data: 

C{pi,P2) = C{K,q) ^ 1 + XeM-RoQo " RWs " RWl), (2-9) 

or, for small q 

C{K, q)^l + X{l- Rlql - R^ql - Rlql). (2.10) 

For a static Gaussian source, these HBT radii {Rq, Rs, Rl) would be independent of \K\, 
and would reveal the spatial extent of the source. For a general dynamic source, the radii 
can be complicated functions of K, and their interpretation more complicated. 
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Chapter 3 
HYDRODYNAMICS 

3.1 Introduction 

Before delving in to the specifics of viscous hydrodynamic models of relativistic heavy ion 
collisions, it is useful to know a bit about the theory of hydrodynamics in general and in 
particular the development of relativistic viscous hydrodynamics (see Paul Romatschke's 
lecture notes |l^ for a nice, more detailed, treatment). 

Hydrodynamics — also known as fluid dynamics — is the theory governing the motion of 
fluids. As the name implies, it was initially developed to describe the dynamics of water, 
but can be applied to fluid behavior of a wide range of materials. It can be thought of as 
an effective theory describing the long wavelength behavior of a system that has sufficient 
separation of scales such that this macroscopic motion is so slowly varying in space and 
time so as to be insensitive to the microscopic dynamics. In the case of water, for example, 
if the macroscopically-averaged quantities such as pressure and temperature change very 
slowly in space compared to the average distance between molecules and very slowly in 
time compared to the scattering rate of the individual molecules, it will behave according 
to the equations of hydrodynamics. 

Likewise, there should be hydrodynamic regimes for systems consisting of a collection 
of hadrons or even a quark gluon plasma. One can imagine a hypothetical large system 
that has had enough time to everywhere come very close to thermal equilibrium, yet still 
has a temperature gradient that slowly varies across the system. This system would behave 
hydrodynamically, even if it consists of a very hot and dense collection of strongly interacting 
matter. 

It is a more difficult question, however, whether the medium created in a relativistic 
heavy ion collision interacts strongly enough to behave like a fluid for the short period 
of time before it flies apart. Although it now appears from the success of hydrodynamical 
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models that this is hkely the case, it was not at all obvious that the hydro dynamic description 
should be correct and it was originally believed by many to be unlikely. 

3.2 Non-Relativistic Fluid Dynamics 

The conventional (non-relativistic) formulation of the hydrodynamic equations describes 
the evolution of the fluid velocity v{t,x), the pressure p(t, x ), a nd the mass density p{t,x) 
of a fluid at each point in space and time via the equations [3], [2^ §2, 



dtp + pd ■ V + V ■ dp = . (3.1) 
dtv + (^v ■ d^ V = —-dp , (3.2) 

These equations are referred to as the continuity equation (|3.ip and Euler equations (|3.2p . 

respectively, and are simply the statements of conservation of mass and momentum for a 

continuous fluid without dissipation (an "ideal" fluid). To close the set of equations another 

relation is needed — usually given as an equation of state of the material p = p{p). The 

Euler equations can be generalized to treat dissipative effects 

9v* ^ ^ 9f * _ 1 dp 1 9n'^* 
dt dx'' p dx^ p dx'' ' 

n^^ = _ + ^ _ ls''^^] -cd^'— . (3.4) 

V dx'' dx'^ 3 dx' J dx' 



2l| 



20(]§15. The coefficients in the viscous 



These are called the Navier-Stokes equations 
stress tensor II'^' are termed the shear viscosity (r/) and bulk viscosity Their values, 
like the form of the equation of state, depend on the specific fluid in question and encode 
information about the microscopic dynamics of that material. 

3.3 Relativistic Ideal Hydrodynamics 

Let us generalize this so that it can be applied to a relativistic system. Any system can 
be characterized by its energy-momentum tensor T^'^{x), which is a symmetric tensor that 
describes the distribution of energy and momentum in the system. In a given reference frame 
the time-time component T^^ is the energy density, the time-space component T*^* = T**^ 
is the i'th component of the momentum density, and the space-space component T*'^ is the 
flux of i'th momentum across the x'' surface. 
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In relativistic notation, the statement of conservation of energy and momentum is simply 

df.T^"' = 0. (3.5) 

Any other conserved quantity (e.g., electric charge, baryon number, etc.) is characterized by 
a conserved current j^{x) that describes its charge density and current. The conservation 
equations then include an additional equation, 

« = 0, (3.6) 

for each conserved quantity (labeled by n). For simplicity it will be assumed that these 
additional conservation equations are unimportant to the motion of the fluid and so will be 
neglected in the following. 

Define the local fluid rest frame at each point in space-time as the zero-momentum 
frame, T^^{x) = 0. The velocity of this local rest frame with respect to a fixed lab frame 
defines a fluid 4-velocity u^{x) such that in the rest frame ^i^est ~ (1; 0)0,0) (recall that for 
a 4-velocity, = u^u^ = 1). So then u^T'^'^ = eu^, where e(x) is defined as the energy 
density in this fluid rest frame. 

The equations of ideal (relativistic) hydrodynamics then emerge from the conservation 
equations 13.51 when the energy-momentum tensor has the property that it is isotropic (i.e., 
rotationally invariant) in the local rest frame: 

/ e \ 

p 
p 

\ p y 

In an arbitrary fixed reference frame, in covariant notation, this is 

TtLi = + P) u'^W' -pg'"' = e u^^u' - p (3.7) 

where p{x) is then the isotropic pressure in the rest frame and g^'^ = diag(l, —1, —1, —1) is 
the metric tensor. A^"^ = [g^^ — u^u'^) is a projection operator on the space orthogonal to 
the fiuid velocity. It has the properties Af'^u^ = A^^^u^, = and A^'^A^ = A''". It is often 



rpt^U 

ideal, rest 
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useful to use this projector to express the hydrodynamic equations projected into directions 
parallel {uudpT^^^) and perpendicular {A"d^T^'^) to the fluid velocity. Explicitly they are: 

^^^^M^dL = + P)^^.u^ + ^^d^^ = (e + P)df^^^ + De = 0, (3.8) 
Kd^TlZ^i = (e + P) u^d^u'' - A^"5^P = (e + p)I)n" - V> = , (3.9) 

where we have also introduced shorthand notation for the projection of derivatives parallel 
{D = u^d^) and perpendicular (V" = A^"5^) to the fluid velocity. 

This system of equations is closed by specifying the equation of state p = p{e). If this 
is the equilibrium equation of state for the system in question, these equations describe 
a system that is everywhere in local thermal equilibrium. Thus, the language used when 
talking about hydrodynamics is often that of thermal equilibrium, but note that all that is 
required for Equations 13.81 and 13.91 to be valid is isotropy in the fluid rest frame. 

For fluid velocities much less than the speed of light [u^ ~ (l?^)]) when the energy 
density is dominated by the mass density (e ~ p), these relativistic ideal hydrodynamic 
equations reduce to the non-relativistic Euler and continuity equations from [section 3.21 
J. 

It should be noted that these (non-linear) ideal hydrodynamic equations contain in- 
stabilities that make them impossible to solve numerically — at least using the most naive 
of algorithms. To avoid these problems, numerical algorithms are used that contain what 

g viscous terms 

to the equations also flxes these instabilities, and so solving the viscous hydrodynamic 
equations will not require these algorithms. 

3.4 Relativistic Viscous Hydrodynamics 

These equations can be generalized to include dissipative, or viscous, effects. This is done 
by allowing for a more general energy- momentum tensor: 



(3.10) 
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U^'^ is the viscous stress tensor that inc 
hydrodynamic equations then become 



udes the contributions to T^'^ from dissipation. The 



1^ 



De + {e + p)d^u^' - n^'^ V(^n,) = , 

(e + p)z)n" - v> + A°a^n^^ = , 



(3.11) 



where the (...) denote symmetrization, e.g. 

1 



(V^-u^ + V 



Of course, one must still specify the form of IT^"^. Determining the correct form to use turns 
out to be less straightforward than in the non-relativistic case, and until recently there has 
been a number of versions in use. 



3.4-1 Relativistic Navier-Stokes equations 

If the macroscopically- averaged quantities from ideal hydrodynamics (e, p, u^) vary slowly 
in space and time, it can be useful to build a controlled gradient expansion of T^'^; i.e., an 
expansion in powers of derivatives of these quantities. Using this perspective, T^^g^^ = Tq"^ 
is just the zeroth order term in such an expansion, and 11'^'^ contains first and higher order 
derivative terms. 

To first order in gradients, there are only two independent terms that are consistent 
with the symmetries of 11^'^ (it must be symmetric, H^'^ = 11'^^^, so that T^"^ also remains 
symmetric, and it must be transverse to the fiuid velocity, u^H^'^ = 0, to retain the definition 
of the zero- momentum rest frame). They are usually separated into a traceless term tt^'^ 
and the remainder 11 

^^,u ^ j,i^u ^ ^f,u ^ ^/.z^jj ^ rjVK"'^ + C A^^'Vau'^. (3.12) 

The angle brackets define a quantity that is traceless, symmetric, and transverse 

2 

or, in general 

A'^^'B^'^ =P°^S^^'B'' , (3.13) 
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with = AqAJ^ + A^A^ — |A'*''Aq(3. Upon plugging into the conservation equations, 
this results in what are termed the relativistic Navier-Stokes equations, since they reduce 
to Navier-Stokes in the non-relativistic limit, ij is then identified as the shear viscosity and 
C the bulk viscosity. Ideal hydrodynamics is recovered when these transport coefficients are 
set to zero (appropriate when gradients are so small compared to the coefficients that the 
terms can be neglected). These terms are often derived in the literature by demanding that 
the .cond law of „a..cs alway. be obeyed, 3,^ > 0. wbe.e = . tbe 

entropy density current [la], rather than using the perspective of a gradient expansion to 
first order. 



3.4-2 Causality restored: second- order relativistic viscous hydrodynamics 

There are problems with the relativistic Navier-Stokes equations, however, which in general 
make them difficult to solve numerically. This is caused by the presence of acausal signal 
propagation, and associated instabilities. 

The problem can be illustrated by considering small perturbations of the energy density 
and fluid velocity, and tracking how these disturbances travel through the medium. If one 
decomposes the perturbations into Fourier modes, one finds that the diffusion speed of a par- 
ticular mode increases linearly with wavenumber. A mode with arbitrarily large wavenumber 
will have an arbitrarily large speed (larger than the speed of light), and causality is violated 

Q. 

On the surface this shouldn't necessarily be worrisome. Large wavenumber (or short 
wavelength) modes are outside the realm of applicability of hydrodynamics — if there is 
significant short distance behavior, the gradient expansion will not converge and hydro- 
dynamics is not an appropriate description of the system anyway. In practice, however, 
this acausal behavior causes instabilities that make constructing a numerical solution with 
arbitrary initial conditions impossible [22| . 

It turns out that, as with the (unrelated) instabilities of ideal hydrodynamics, these 
instabilities can be removed by adding higher order gradient terms to the equations. 

At second order in gradients, one can construct 15 independent terms (in addition to 
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the zeroth and first order terms already mentioned) [231]. In the shear (traceless) sector we 
have 



K 

+ 2 



- K*UaUi3R''^^"'^^ - r?r;^(V • u)a^'^ + ^V^^lnsV''^ Ins , (3.14) 

while in the bulk sector the most general form is 

n = C (V • n) - CmD (V-u)- iia^^'a^, - 6 (V • uf 

- iz^^'^'^i.v + e4V^lnsV^lns + ^5^? - i^u^vPR^fi , (3.15) 

where s is the entropy density. Since we are assuming vanishing conserved charge densities 
(i.e., zero chemical potential), this is given by s = Here we have also introduced the 
notation o^^ = V'^^u'^\ and the fluid vorticity is defined as uj^^ = — V'^'n'^]. Square brackets 
indicate an antisymmetrized quantity 

In a general curved space-time the Riemann tensor is non-zero 

oA _ Q pA Q -pA I -pre pA pre pA 

with r^j^ = \g^^ {d^gpi, + di^gpf^ — dpg^^) , as well as the Ricci tensor Rp^ = R^^^u ^^'^ Ricci 
scalar R = R^. 

Therefore, without any additional information, there are in principle 15 possible second- 
order transport coefficients multiplying these terms: Tt^, r*, k, k*, Ai, A2, A3, A4, rn, ^1, ^2, 
Csj ^4) ^5 aiid ^61 in addition to the first order transport coefficients rj and (. If one is lucky, 
these transport coefficients can be calculated from the underlying microscopic theory (e.g., 
QCD in the case of relativistic heavy ion collisions). Often, however, this is not feasible (as 
is the case for QCD at the moment, although there is ongoing work). If they cannot be 
computed from first principles, the transport coefficients can be treated as free parameters 
that are then constrained by experimental data. In this case, it becomes important to 
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understand which terms are necessary (or appropriate) to keep in a hydrodynamic simulation 
of a given problem. 

One can use general arguments to reduce the number of free parameters: requiring the 
positivity of the divergence of the entropy current provides 2 extra constraints, leaving 
only 13 completely independent transport coefficients [23] . Also, when space-time can be 



assumed flat, the terms involving the Riemann tensor and the Ricci tensor and scalar 
drop out, reducing the number of independent second order transport coefficients to 11, in 
addition to the 2 first order coefficients. 

For most problems, however, this situation is still not satisfactory. Having so many free 
parameters (as well as the fact that their temperature dependence is a priori unknown) 
would significantly reduce the predictive power of any simulation. On the other hand, while 
it is true that one can eliminate the problems of the relativistic Navier-Stokes equations by 
adding only a single second order term, the precise form of the term used can affect how 
sensitive the results are to the value of the corresponding transport coefficient (as well as 
affecting the interpretation of the resulting value in the context of the underlying theory), 
and so one must be convinced that keeping only that part of the second order expansion is 
justified. 



3.4-3 Miiller- Israel- Stewart theory 



This issue of a causal set of relativistic viscous hydrodynamic equations was originally 
studied in depth by Miiller 
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271, 



28|. Thus, 



25[ , and separately by Israel and Stewart 
many of the second order extensions to viscous hydrodynamics are commonly referred to as 
Miiller-Israel-Stewart theory, although the exact form of the equations that have been used 
varies. An extensive history of the various forms used and derivations thereof is beyond the 
aim of this summary, although some generic comments can be made. 

In all cases one has a term whose coefficient acts as a relaxation time which, when 
it is larger than a certain value, eliminates any possible acausality. In the shear sector 
this relaxation time is usually called r^r, and it typically multiplies some combination that 
includes the term ^Da^'^^ (e.g. in Equation 3.14^ . 
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Besides causality, another important consideration is the second law of thermodynamics, 

(3.16) 



d^s^" > . 



which should be obeyed in any physical system. In the absence of dissipation, the entropy 
current is given by = su^. In ideal hydrodynamics, the inequality is saturated, and 
the entropy does not increase with time. Assuming this form for the entropy current, one 



can show that the relativistic Navier-Stokes equations always obey Equation 3.16 for any 
(non-negative) value of the transport coefficients jl8l ]. 



When there is dissipation, however, the entropy current can also have gradient terms. 
If one assumes that the entropy current has to be algebraic in the hydrodynamic degrees 
of freedom, and that the deviations from equilibrium are small enough that higher order 
corrections can be neglected, it can be shown that the entropy current has to be of the form 



25 



(3.17) 



One can then plug this into Equation 3.16| and derive a form for the second order hydro- 
dynamic terms that is then guaranteed to always obey the second law of thermodynamics 



18l |. There is doubt, however, that the assumption that the entropy current has to be 



algebraic in the hydrodynamic degrees of freedom is necessarily true and therefore also that 



Equation 3.17 is the really the most general form 
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301 ]. Also, it may be an unnecessarily 



strong requirement to demand that the structure of the hydrodynamic equations must be 



such that Equation 3.16 is satisfied in any regime and for any combination of values of 
the transport coefficients — it is only necessary that the second law of thermodynamics be 
obeyed in actually physically realized (or at least realizable) systems. 

Another tack that is often taken is to look to the kinetic theory of gases for guidance. 
Kinetic theory is a description of a system in terms of collisions of dilute particles (or quasi- 
particle states that are long lived compared to the scattering rate). By considering small 
departures from equilibrium, one can derive another version of Miiller-Israel-Stewart theory 



la ]. Unfortunately, although kinetic theory has a large range of applicability and can thus 



often provide much insight into various physical behavior, such a particle-based description 
is not valid for, e.g., the strongly coupled non-abelian plasma which may exist for part of 
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the evolution of a relativistic heavy ion collision system. Therefore, it is uncertain that the 
resulting second-order viscous hydrodynamic equations are general enough to describe the 
evolution of such a system. 

3.4-4 Conformal relativistic viscous hydrodynamics 

If there is reason to believe that a system is approximately scale- invariant, the form of the 
hydrodynamic equations are greatly restricted. To be precise, if one assumes a conformal 
symmetry in the underlying physics, the total number of possible independent first and 
second order transport coefficients is reduced to 6 — and only 4 in the case of flat space. 

The conformal group is the set of symmetry transformations that consists of scale trans- 
formations and special conformal transformations. A theory that is conformally symmetric 
(a.k.a. "conformal" ) has the property that its action is invariant under a Weyl transforma- 
tions of the metric, 



where w{x) is an arbitrary function of x. 

This implies that (to second order in derivatives) the trace of the energy-momentum 
tensor vanishes, and that T^'^ transforms under a Weyl rescaling as 



9tiu 



9tiu = e 



2w{x) 



(3.18) 




,6ui(x) 



(3.19) 



Imposing these conditions gives as the most general form 311] 





(3.20) 



Note in particular that there is no bulk viscosity. To this order in derivatives, it is valid to 
replace rja^^^ — > tt'^'', and so the form that is more convenient to actually solve numerically 



22 



IS 



K 

+ 2 
Ai 



A. 



2ry^ 2?? 2 



(3.21) 



In the expectation that the medium created in a relativistic heavy ion colhsion is approxi- 
mately conformal, this is the form that will be used in the viscous hydrodynamic simulations 
of the following chapters. 
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Chapter 4 

MODELING HEAVY ION COLLISIONS USING VISCOUS 

HYDRODYNAMICS 

4.1 Anatomy of a Heavy Ion Collision 




Figure 4.1: Space-time cartoon diagram of a heavy ion colhsion, indicating the pre- 
equihbrium stage, hydrodynamic stage, and the system after it has frozen out. Here it 
is assumed that there is hydrodynamic evolution on both sides of the phase transition as 
weh as approximately boost invariant evolution. 

Even if a heavy ion collision system does behave hydrodynamically for a significant pe- 
riod of its evolution, that is not the whole story, of course. Once the proper hydrodynamic 
equations are set, one must still specify the boundary conditions. There is a finite period 
of time at the beginning of a collision before which the system can equilibrate (or at least 
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isotropize — recall Isection 3.3P and begin behaving hydro dynamically. Likewise, as the sys- 
tem expands and cools, there will be some point at which the system no longer interacts 
strongly enough for hydrodynamics to be a valid description. Eventually, the particles get so 
far apart that they completely cease interacting and these free particles are what ultimately 
get detected. Therefore one must define the process by which the system "freezes out". 
Once the initial conditions and freeze out algorithm have been specified, the simulations 
can be run numerically and experimental observables calculated. Details of these choices 
will be given later in this chapter when describing the viscous hydrodynamic simulations 
done by the author in collaboration with Paul Romatschke. 

Ideal hydrodynamic simulations were previously done, however, offering a remarkably 



good description of the experimental data for bulk properti es ( multip^ 



3, 
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i city . 



34], 



radial and elliptic 
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flow) of low pt particles for heavy-ion collisions at RHIC 

Upon closer inspection, however, not all of this success can be attributed to modeling 
the system as an ideal fluid. For instance, the energy density distribution which is used 
as an initial condition for the hydrodynamic equations is customarily chosen such that the 
output from the hydrodynamic model matches the experimental data for the multiplicity. 
Furthermore, the time where the hydrodynamic model is initialized as well as the tem- 
perature (or energy density) at which the hydrodynamic evolution is stopped are typically 
chosen such that the model output matches the experimental data for the radial flow. After 
these parameters have been flxed, only the good description of experimental data for the 
elliptic flow coefflcient can be considered a success for ideal hydrodynamics (in the sense 
that it is parameter- free). 

In order to make progress and learn more about the properties of matter created at 
RHIC, the task is now to both test and improve this ideal hydrodynamic model. The obvious 
framework for this task is dissipative hydrodynamics, since it contains ideal hydrodynamics 
as the special case when all dissipative transport coefficients (such as shear and bulk viscosity 
and heat conductivity) are sent to zero. If the value of the transport coefficients were known 
(e.g. by some first principle calculation), then one could use dissipative hydrodynamics 
to constrain e.g. the initial energy density distribution, which is chosen conveniently in 
the ideal hydrodynamic models. Or otherwise, choosing again physically acceptable initial 
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conditions, one is able to constrain the allowed ranges of the transport coefficients. Despite 

" 

recent progress in first principle calculations 
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37 



38|, 



3S 



401, 



41 



42 



43|, |44|, |4|, the 



values of the hydrodynamic transport coefficients for QCD in the relevant energy range are 
poorly constrained to date, so the second option is currently the only viable possibility. 



For RHIC, the first step in this direction was carried out by Teaney [46(], who provided 
estimates for the sign and size of corrections due to shear viscosity. This famous calculation, 
however, did not provide a description of experimental data for non-zero viscosity, because 
it was not dynamic and the initial conditions could not be altered. Only very recently, the 
first hydrodynamic calculations with shear viscosity describin g p article spectra for central 



and non-central collisions at RHIC have became available 
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49j. 



Several other 
matching to data 



groups 



50, 



51 



ha ve p roduced numerical codes capable of performing similar 



52, 



53|, 



54|, 



However, the precise formulation of the viscous hydrodynamic equations themselves has 



long been debated (recall Isection 3.4p . For the case of non- vanishing shear viscosity only, it 



was shown recently 



3ll | that the most general form implies five independent terms of second 



order in gradients. This form is general enough to describe the hydrodynamic properties 
of (conformal) plasmas both for weakly coupled systems describable by the Boltzmann 
equation as well as infinitely strongly coupled plasmas, which are accessible via Maldacena's 
conjecture [58 1. 



The aim of this chapter is to now apply this new set of equations for relativistic shear 
viscous hydrodynamics to the problem of heavy-ion collisions at RHIC. In Isection 4.2) we 
review the setup of conformal relativistic viscous hydrodynamics and our numerics for the 



simulation of heavy-ion collisions. In lsection 4.31 details about the two main models of initial 
conditions for hydrodynamics are given. Section [4.41 contains our results for the multiplicity, 
radial flow and elliptic flow in Au-|-Au collisions at top RHIC energies, as well as a note on 



the notion of "early thermalization" . We conclude in Isection 4.51 
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4.2 Setup 

We use the most general form of the second order viscous hydrodynamic equations, which 
as a reminder are given by (see [section 3.4p 



K 

+ 2L- 



2r]^ 2r] 2 



(4.1) 



The coefficients t^, k, Ai, A2, A3 are the five new coefficients controlling the size of the allowed 
terms of second order in gradients. Having an application to the problem of heavy- ion 
collisions in mind, the above set of equations can be simplified: for all practical purposes 
spacetime can be considered flat, such that both the Riemann and Ricci tensors vanish 
identically. Thus, only the four coefficients r^, Ai, A2, A3 enter the problem. 



4.2.1 A note on bulk viscosity and conformality 

Besides shear viscosity, QCD also has non-vanishing bulk viscosity (" which can be related 
to the QCD trace anomaly [sO] 

C^Tli = e- 3p. (4.2) 

QCD lattice simulations seem to indicate that the ratio bulk viscosity over entropy density 
■s, C/s, is small compared to r]/s except for a small region around the QCD deconfinement 



60 



61 



62l |. If we are interested in 



transition temperature, where it is sharply peaked 
describing effects from shear viscosity only, we are led to consider C = 0, or conformal fluids, 
This has been the main guiding principle in Ref. 



311 ] and as a consequence Equation 4.1 
obeys conformal invariance, unlike most other second-order theoriej^. 



4.2.2 First steps: 0+1 dimensions 

In order to get a crude estimate of the effect of viscous corrections, let us consider the 
arguably simplest model of a heavy-ion collision: a system expanding in a boost-invariant 



^Note that Muronga derived a version of |Equation 4.1| in Ref. [29y that turns out to obey conformal 
symmetry. 
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fashion along the longitudinal direction and having uniform energy density in the transverse 



plane. Introducing the Milne variables proper time r = yji^ — and space-time rapidity 
^ = arctanh(2:/t), boost invariance simply translates to requiring all hydrodynamic vari- 
ables (e, u'^, n'^'^) to be independent of rapidity, and tensor components n^,n^^ to vanish. 
Assuming uniformity in the transverse plane furthermore requires independence from the 
transverse coordinates = {x,y). Even though this means that all the velocity compo- 
nents except u'^ are zero, the system is nevertheless non-trivial in the sense that the sum 
over velocity gradients does not vanish, V ^u^ = ^, sometimes referred to as "Bjorken flow". 

In a way one has modeled an expanding system in static space-time by a system at 
rest in an expanding space-time. This has been achieved by transforming to the Milne 
coordinates r, ^, where the metric is g^jy = diag{grT , dxx , 9yy , 9^^) = (l, — l, — l, —t"^)- Note 
that even though the spacetime in these coordinates is expanding, it is nevertheless flat (e.g. 
has vanishing Riemann tensor). 

In this 0+1 dimensional toy model, the viscous hydrodynamic equations become excep- 
tionally simple 3l|, 



8.. = -i±i! + 5t 

r r 

a^ul = -5l + J!L_^^|--^(^|)^ («) 

The Navier-Stokes equations are recovered formally in the limit where all second-order 
coefficients vanish (e.g. t^, Ai — > 0); then, one simply has 

n« = ^. (4.4) 



The equations (j4.3p can be solved numerically along the lines of 63|,|6j]. At very early times, 
where n| > (e + p), the Navier-Stokes equations indicate an increase in energy density and 
a negative effective longitudinal pressure p — n|. Since gradients V^n'^ = l/r are strongest 
at early times, this suggests that one is applying the Navier-Stokes equations outside their 
regime of validity. Theories including second order gradients may be better behaved at early 
times, but eventually also have to break down when gradients become too strong. Here we 
want to study the effects of the second order coefficients on the value of the shear tensor at 
late times, where a hydrodynamic approach should be valid. 
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To this end, let us study the deviation of the shear tensor from its first order value, 
611 = rig — 1^. At late times. Equation 4.3 implies e ~ t~^/^, so ~ r^^. Thus, if 511 is 



small compared to the first order value, from Equation 4.3 we find 

6U 



4?? /2t^ 2Ai\ 
3r \ 3r Stt] J 



(4.5) 



For a strongly coupled = 4 plasma 
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65|, 



66], one hai 



1 
s 



1 

4^' 



In 2 



Ai 



V 



(4.6) 



2ttT ' 27rT' 

and thus n| is larger than its first order value by a factor of 1 + g"!^^^ . For RHIC, Tr > 1 
is a reasonable estimate, so one finds that the second order corrections to n| increase its 
value by a few percent over the first order result. 

As an example on the importance of obeying conformal invariance, imagine dropping 
the term involving Vqu" in the first line of [Equation 4.1 Redoing the above calculation 
one finds 

(4.7) 



6t \ T Stt] J 



which indicates a nearly ten-fold increase of the size of 611 for the non-conformal theory. 
For a weakly coupled plasma well described by the Boltzmann equation jsi], where one has 



6r; 



(Ai is unknown but generally set to zero in Miiller- Israel-Stewart theory ), the 



effect may be less pronounced, but still one qualitatively expects second-order effects to be 
anomalously large if conformal invariance is broken in an "ad- hoc" manner. 

Clearly, the above estimates are not meant to be quantitative. Indeed, even the sign 
of the correction may change when allowing more complicated (e.g., three-dimensional) 
dynamics. However, the lesson to be learned from this exercise is that second-order gradients 
can and indeed do modify the shear tensor from its first order (Navier-Stokes) value. This 
is physically acceptable, as long as the second-order corrections are small compared to 
the first order ones (otherwise the system is probably too far from equilibrium for even a 
hydrodynamic description correct to second order in gradients to be valid). A practical 
means for testing this is calculating physical observables for different values of the second- 



^For completeness, we also mention the results k = A2 = — ^ , A3 = from [3ll.[65 
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order coefficients and making sure that the results do not strongly depend on the choice for 
these specific values. 



4-2.3 Including radial flow: lessons from 1+1 dimensions 

Some more insight on the effect of viscous corrections may be gained by improving the 
model of the previous subsection to allow for radially symmetric dynamics in the transverse 
plane (but still assuming boost invariance). This is most easily implemented by changing 
to polar coordinates {x,y) — > (r, 0) with r = y^x^ + and = arctan(y/2;). In this case, 
the only non- vanishing velocity components are and u'\ and hence the vorticity lo^^'^ 
vanishes identically. Although non-trivial, the radially symmetric flow case is still a major 
simplification over the general form Equation 4.1, since again the terms involving k,A2,A3 



drop out. 

Such a formulation allows both for important code tests 16711 as well as realistic sim- 

n 



ulations of central heavy- ion collisions [47[ (note that truncated versions of Equation 4.1 
were used in these works). The advantage of this formulation is that since the equations 
are comparatively simple, it is rather straightforward to implement them numerically and 
they are not very time consuming to solve since only one dimensional (radial) dynamics 
is involved. The shortcoming of simulations with radially symmetric flow profiles ("radial 
flow") is that by construction they cannot be matched to experimental data on the impact- 
parameter dependence of multiplicity. Thus, the considerable freedom in the initial/final 
conditions inherent to all hydrodynamic approaches cannot be eliminated in this case. 

For this reason, we will choose not to discuss the case of radial flow here in more detail, 
but rather will comment on it later as a special case of the more general situation. 



4--2.4 Elliptic flow: 2+1 dimensional dynamics 

Retaining the assumption of boost invariance, but allowing for general dynamics in the 
transverse plane, it is useful to keep Cartesian coordinates in the transverse plane, and 
thus n'^,u^',u^ are the non- vanishing fluid velocities. The main reason is that e.g. in polar 
coordinates the equations for the three independent components of 11'^'' would involve some 
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extra non- vanishing Christoffel symbols (other than T^^ = t,T^^ = l/r). 



Fortunately, the case of two dimensions is special insofar as the only nontrivial compo- 
nent of the vorticity tensor, namely to^y, fulfills the equation 4S| 



■ UJ 



xy 



Dp Du" 



e + p 



o(n3). 



(4.8) 



which can be derived by forming the combination V^Du^ — V^Du^. The expression C(n^) 
denotes that the r.h.s. of Equation 4.8| is of third order in gradients, and thus should 
be suppressed in the domain of applicability of hydrodynamics. For heavy-ion collisions, 
typically ^u^ > ^, so that for an equation of state with a speed of sound squared = 



^^^^ ~ ^, [Equation 4.8 translates to ^^xy < unless Dlnn"^ > (1 — c^)V^n^'. In particular, 
this implies that in general uj^y = is a stable fix-point of the above equation and hence 
we expect uj^^ to remain small throughout the entire viscous hydrodynamic evolution if it 
is small initially. 



Generically, one uses u 



x,y 



as an initial condition for hydrodynamics 



681 ]. which 



implies = initially. Therefore, to very good approximation we can neglect the terms 



involving vorticity in Equation 4.1, such that again only the second-order coefficients t^^, Ai 
have to be specified. 



The equations to be solved for 2-1-1 dimensional relativistic viscous hydrodynamics are 
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then (in components) 



3U'^ TjrU'^ TjrU^ 



u 



2 

3^ 

V^gu^) = 2r^A««rJ^n^ - ^r^A^^VX- (4-9) 



Here and in the following Latin indices collectively denote the transverse coordinates x, y and 
the relation n^II'^'' = has been used to derive the above equations (similarly, u^V (^^Uy\^ = 
can be used to obtain the other non-trivial components needed). Note that this particular 



form of Equation 4.9 has not been simplified further since it roughly corresponds to the 



equations implemented for the numerics of 



4a], and is meant to facilitate understanding of 



the code |69i]. A simple algorithm to solve [Equation -£9 has been outlined in [67] and will 



be reviewed in the next subsection for completeness. 

4--2.5 A numerical algorithm to solve relativistic viscous hydrodynamics 

The first step of the algorithm consists of choosing the independent degrees of freedom. For 
boost-invariant 2+1 dimensional dynamics, a sensible choice for this set is e.g. e, u^, u^, H^^, 
U^y , . The pressure is then obtained via the equation of state p(e), and the only other 



non-vanishing velocity as n"^ = ^1 + n| + u^. Similarly, the other nonzero components of 
n^'^ are calculated using the equations 11^ = 0, u^H^'^ = 0. 

Given the value of the set of independent components at some time r = tq, the aim 
is then to construct an algorithm from [Equation 4.9 such that the new values of the set 
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can be calculated as time progresses. Note that in Equation 4.9, time derivatives of the 



independent component set enter only linearly. Therefore, Equation 4.9 may be written as 
a matrix equation for the derivatives of the independent component set, 

/ ann n.m rtriK \ / (^7-6 A / hn \ 

drU^ 



«00 ^01 
«10 Oil 



«05 

ai5 



bo 
bi 



(4.10) 



y aso 051 ... 055 J \ dr^y J \ be / 
Denoting the above matrix and vector as A and b, respectively, a straightforward way to 
obtain the time derivatives is via numerical matrix inversion, 



\ drUyy J 

Choosing a naive discretization of derivatives 



• b. 



(4.11) 



dxf{x) 



f{x + a) - f{x - a) 



(4.12) 



6t ' 2a 

which is first order accurate in the temporal grid spacing 5t and second order accurate in 
the spatial grid spacing a, one can then directly calculate the new values of the independent 



component set from Equation 4.11[ 

Note that for ideal hydrodynamics, the algorithm [Equation 4.1l| would fail for this naive 



discretization 



70( 1 . The reason is that ideal hydrodynamics is inherently unstable to high 



wavenumber fluctuations (which can be thought of as the basis for turbulence). For ideal 
hydrodynamics, one thus has to use a discretization which amounts to the introduction of 
numerical viscosity to dampen these fluctuations. Luckily, viscous hydrodynamics does not 
suffer from this problem because it has real, physical viscosity inbuilt. It is because of this 



reason that the naive discretization can be used in the algorithm Equation 4.11 without 
encountering the same problems as in ideal hydrodynamics, as long as a finite value for the 
viscosity rj is usecj^. While applicable to sufficiently smooth initial conditions, the above 



^ In practice we have used 2 > 10 Typically, between ^ = 10 ^ and ^ — 10 there are no significant 
changes in the hydrodynamic results and we refer to = 10"'' as "ideal hydrodynamics" . 
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algorithm is too simple to treat strong gradients such as the propagation of shocks, and 
should be replaced by a more involved scheme in such cases. 

Since matrix inversions are computationally intensive, one can speed up the numerics 
by expressing dr^^^'^ in terms of drU^, dr^- Inserting these in the equations for Du^ and De, 
one only needs to invert a 3 x 3 matrix to obtain the new values of the energy density and 



fluid velocities. This approach has been used in 



^,1481,16 



4-2.6 Initial conditions and equation of state 

As outlined at the beginning of the chapter, any hydrodynamic description of RHIC physics 
relies on given initial energy density distributions. Two main classes of models for boost- 
invariant setups exist: the Glauber models and the Color-Glass-Condensate (CGC) models. 

As will be shown in the following, both model classes can give a reasonable description 
of the experimentally found multiplicity distribution, but they differ by their initial spatial 
eccentricity. A detailed discussion of the initial conditions will be given in subsequent 
sections. 

Besides an initial condition for the energy density, one also needs to specify an initial 
condition for the independent components of the fluid velocities and the shear tensor. For 
the fluid velocities we will follow the standard assumption that these vanish initially [s^. 
Finally, when using the set of equations (j4.9p . one also has to provide initial values for the 
independent components of U'^'^. Extreme choices are H'^'^ = and a shear tensor so large 
that a diagonal component of the energy-momentum tensor vanishes in the local rest frame 
(e.g. = p, or zero longitudinal effective pressure), with the physical result expected 



somewhere in between (see e.g. the discussion in j7ll|). 

Once the initial conditions for the independent hydrodynamic variables have been speci- 
fied, one needs the equation of state to solve the hydrodynamic equations (j4.9p . Aiming for a 
description of deconfined nuclear matter at zero chemical potential, a semi-realistic equation 
of state has to incorporate evidence from lattice QCD calculations [3] that the transition 
from hadronic to deconfined quark matter is probably an analytic crossover, not a first or 
second order phase transition as often used in ideal hydrodynamic simulations. On the 
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other hand, continuum extrapolations for the value of the energy density and pressure for 
physical quark masses are still not accessible with high precision using current lattice meth- 
ods. For this reason, we will employ the equation of state by Laine and Schroder [3], which 
is derived from a hadron resonance gas at low temperatures, a high-order weak-coupling 
perturbative QCD calculation at high temperatures, and an analytic crossover regime inter- 
polating between the high and low temperature regime, respectively. For hydrodynamics, 
an important quantity is the speed of sound squared extracted from the equation of state, 



^2 ^ For completeness, we reproduce a plot of this quantity in Figure 4.2 



^.2.1 Freeze- out 

At some stage in the evolution of the matter produced in a heavy-ion collision, the system 
will become too dilute for a hydrodynamic description to be applicable. This "freeze-out" 
process is most probably happening gradually, but difficult to model realistically. A widely- 
used approximation is therefore to assume instantaneous freeze-out whenever a certain fluid 
cell cools below a certain predefined temperature or energy density (see Q, 13] for differ- 
ent approaches). The standard prescription for this freeze-out process is the Cooper- Frye 
formula [t^, which allows conversion of the hydrodynamic variables (energy density, fluid 
velocity,...) into particle distributions. 
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Specifically, in the case of isothermal freeze-out at a temperature Tf, the conversion 
from hydrodynamic to particle degrees of freedom will have to take place on a three- 
dimensional freeze-out hypersurface S, which can be characterized by its normal four-vector, 

The spectrum for a single particle 



and parametrized by three space-time variables H] ■ 
on mass shell with four momentum p^^ = (E, p) and degeneracy d is then given by 

d^N d 



E- 



(4.13) 



d^p (2^)3 

where dT,^ is the normal vector on the hypersurface S and / is the off-equilibrium distri- 
bution function. 

Originally, the Cooper-Frye prescription was derived for systems in thermal equilibrium, 
where / is built out of a Bose or Fermi distribution function /o (x) = (exp[(x) ± l]-i), 
depending on the statistics of the particle under consideration. In order to generalize it to 
systems out of equilibrium, one customarily relies on the ansatz used in the derivation of 
viscous hydrodynamics from kinetic theory 7a], 



f{x^,pn = fo 



+ Jo 



iT/o 



2T^e + p)' 



(4.14) 



(4.15) 



T J \ T 

For simplicity, in the following we approximate foix) ~ exp(— x), so similarly 

/ (x^p'^) = exp i-p.u'^/T) [l + 1^^^ • 

The effect of this approximation will be commented on in the following sections. 

In practice, for boost-invariant 2+1 dimensional hydrodynamics, the freeze-out hyper- 
surface = (S*, S^, S^) = (t, X, y, z) can be parametrized either by r, ^ and the polar 
angle (j), or by x, y, ^: 



t = T cosh ^ 
X = x{t, (j)) 

y = yi-r, (p) 

z = T sinh ^ 
The normal vector on is calculated by 



t = r(x, y) cosh ^ 

X = X 

y = y 

z = r(x, y) sinh^ 



(4.16) 



as^ 



dr 



drdcfxi^ 



(( dx dy dy dx \ dy dx ( dx dy dy dx \ \ 
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Figure 4.3: Space-time cut through the three-dimensional hypersurface for a central collision 
within the Glauber model. Simulation parameters used were a = 1 GeV~^, tq = 1 fm/c, 
Tj = 0.36 GeV, Tj = 0.15 GeV, = 6^ and Ai = (see next sections for definitions). 
As can be seen from the figure, inclusion of viscosity only slightly changes the form of the 
surface. 



and similarly for the other parametrization {t^ . 

For a realistic equation of state, at early times the freeze-out hypersurface will contain 
the same transverse coordinate values {x, y) for different times r (see |Figure 4.3 ). Therefore, 



the parametrization in terms of (x, y, ^) cannot be used for early times. On the other hand, 
the parametrization in terms of (r, 0, ^) contains derivatives of (x, y) with respect to r, which 
become very large at late times (see [Figure 4. 3D . Numerically, it is therefore not advisable 
to use this parametrization at late times. As a consequence, we use the one parametrization 
at early times but switch to the other parametrization at late times, such that the integral 
in [Equation 4.13] is always defined and numerically well-behavecj^. 



In order to evaluate the integral (|4.13p . it is useful to express also in Milne coordinates, 

pfi ^ [p^p^^pV^pi) = (mTCosh(y-^),p^,p^ — sinh(y-e)), (4.17) 

r 



may be possible that other parametrizations may turn out to be more convenient. For instance, it 
is conceivable that performing a triangulation of the three-dimensional hypersurface and replacing the 
integral in (|4.13|l by a sum over triangles could turn out to be numerically superior to our method. 
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where itit = yjm? + p"^ = ~ Pz- Here and in the following Y = aicta.nh.{p^ / E) 

is the rapidity, and m is the rest mass of the particle under consideration. Then the ^ 
integration can be carried out analytically using 



- / cosh'' {Y - OeKp[-zcosh{Y - C)] = i-lTd^Ko{z) = K{n,z), 

2 J-OO 



(4.18) 



where Kq{z) is a modified Bessel function. One finds 
d^N 2d 



d^ 



P 



dTd(l) exp [(^j^u^ + pVuV) /T] x 



( l^l! - 1^1? ) {T^K{l,mTU^/T)+T2K{2,mTU^/T) + T^K{'i,mTU-/T)) 



dr I 



T2 



P'^^-P^'in) {TiK{0,mTuyT)+T2K{l,mTuyT)+T3K{2,mTuyT)) 



1 + 



-2mT 



2T'^{e + p) 
2r2(e + p) ' 



for the (r, i;^), ^) parametrization, and a similar result for the other parametrization of the 
hypersurface. The remaining integrals for the particle spectrum have to be carried out 
numerically unless one is considering the case of a central collision [s^] where the integral 
has an additional symmetry in (p. 

For the simulation of a heavy-ion collision, one then also needs to take into account the 
feed-down process of particle resonances that decay into lighter, stable particles 



80, 



81|. 



Therefore, we calculate the spectra for particle resonances with masses up to ~ 2 GeV and 
then use available routines from the AZHYDRO package [s^ to determine the spectra of 
stable particles including these feed-down contributions. Ultimately, one would be interested 
in describing the last stage of the evolution by coupling the hydrodynamics to a hadronic 



cascade code 



83 



84 



85, 



861 ]. We leave this for future work. 



The particle spectra E 



including feed-down contributions can then be used to 



calculate experimental observables at central rapidity y = , such as radial and elliptic 
flow coefficients, fo,i'2, respectively, and the mean transverse momentum (pt), as defined 



in [section 2.2.11 
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4-2.8 Code tests 



It is imperative to subject the numerical implementation of the relativistic viscous hydrody- 
namic model to several tests. The minimal requirement is that the code is stable for a range 
of simulated volumes and grid spacings o, such that an extrapolation to the continuum may 
be attempted (keeping the simulated volume fixed but sending a — > 0). Our code fulfills 
this property. 



Furthermore, one has to test whether this continuum extrapolation corresponds to the 
correct physical result in simple test cases. One such test case is provided by the 0+1 
dimensional model discussed in [section 4.2.2[ Using initial conditions of uniform energy 
density in the 2+1 dimensional numerical code, the temperature evolution should match 
that of [Equation 4.3 for which it is straightforward to write an independent numerical 
solver. Our 2+1 dimensional code passes this test, for small and large rj/s and different 
values for t-j^, Xi. 



The above test is non-trivial in the sense that it allows to check the implementation 
of nonlinearities in the hydrodynamic model. However, it does not probe the dynamics of 
the model, since, e.g., all velocities are vanishing. Therefore, another test that one can 
(and should!) conduct is to study the dynamics of the model against that of linearized 
hydrodynamics (this test was first outlined in Ref. [67[; see [871 ] for similar considerations). 
More specifically, let us consider a viscous background "solution" with = but non- 
vanishing e(r), n|(T) obeying |Equation 4.3 To first order in small fluctuations 6e, 5u^ , 5Ii^^ 
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around this background the set of equations (14. 9p become 



c%e + ^drUl + An| + (e + p + ln«)a. 



(5u^ + c^S^iJe + didW" 



+ c^a^Je + diSW' = 



9,+ 



5e + 



(e + p) + -n| 



T ^ 







"4 1 




_3r 


(5n^ - 



3r 



+ ^ — n^"" 



3r^r 4r^ 



3tt^ 4r^ « 

— + —djw^ + 
e r,r 



+ 



2r/ 



3r 



3r^ 3 € 

+ In"- 
3r7r 3 

7? 



ai(5n' = 
a-^n^ = 



0, 



(4.20) 



where we have put Ai = and assumed a constant for simphcity. Noting that Sliyy = 
5n| - 5W from ^H^l = 0, |Equat ion 4.20| are a closed set of hnear, but coupled differential 



equations for the fluctuations 6e,6u^ ,6uy ,SIlf,6Il^^ ,611^^ . Doing a Fourier transform, 



6e{T,x,y) 



(4.21) 



(and likewise for the other fluctuations), [Equation 4.20| constitute coupled ordinary differ- 
ential equations for each mode doublet k = (A;^,/c^), which again are straightforward to 



solve with standard numerical methods 

A useful test observable is the correlation function 

{6e{T, xi)5e(T, X2)) 



69( 1 (and analytically for ideal hydrodynamics). 



/(r, xi,X2) 



(4.22) 



where () denotes an ensemble average over initial conditions 6e\^^^^. In particular, let us 
study initial conditions where 6e is given by Gaussian random noise with standard deviation 
A, 

/(ro,xi,X2) = A252(^^_^^) (423) 

and all other fluctuations vanish initially. These initial conditions are readily implemented 
both for the full 2+1 dimensional hydrodynamic code as well as for the linearized system 



Equation 4.20 As the system evolves to finite time r, both approaches have to give the 
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k [GeV] 



Figure 4.4: The correlation function /(r, k) as a function of momentum k = |k| for a lattice 
with a = 1 GeV~^, 64^ sites and averaged over 30 initial configurations (symbols), compared 
to the result from the linearized hydrodynamic equations (lines). 



same correlation function / as long as the linearized treatment is applicable, and hence 



Equation 4.20 can be used to test the dynamics of the full numerical code. 



In practice, note that for the above construction / can only depend on the difference of 
coordinates, 

(,5e(T,xi),5e(r,X2)) _ , _ f 

and therefore in Fourier-space 



/(r,x,-X2) = j ALe"<(«-«)/(r,k) (4.24) 



In the full 2-|-l dimensional numerical code which is discretized on a space-time lattice, 
(5^(k') is regular for any finite a, and one can maximize the signal for /(r, k) by calculating 



the r.h.s. of Equation 4.25 for k' = 0. Similarly, one solution 5e(r, k) per k mode is sufficient 
calculate /(r, k) for the linearized system [Equation 4.20 



The above initial conditions imply /(r = ro,k) = A^, but for finite times characteristic 
peaks develop as a function of |k|, whose position, height and width are sensitive to the 
values of c^, Ttt, r//s and of course the correct implementation of the hydrodynamic equations. 



41 



The comparison between full numerics and linearized treatment shown in Figure 4.4 suggests 
that our code also passes this tesjf]. 

Finally, for the case of ideal hydrodynamics, analytic solutions to the hydrodynamic 



equations are known 



8S 



901 1 ■ Specifically, the code for central collisions 



found to agree with the results from Ref. 



67l | has been 



for ideal hydrodynamics. Since our code agrees 
with Ref. j67|] in the case of central collisions and when dropping the appropriate terms in 
the equations (j4.ip . this provides yet another test on our numerics. 

To summarize, after conducting the above tests we are reasonably confident that our nu- 
merical 2+1 dimensional code solves the relativistic viscous hydrodynamic equations (14. 9p 
correctly. This completes the setup of a viscous hydrodynamic description of relativistic 
heavy-ion collisions. In the following sections, we will review comparisons of viscous hydro- 
dynamic simulations to experimental data, for both Glauber and CGC initial conditions. 

4.3 Initial Conditions: Glauber Model vs. CGC 



4-3.1 The Glauber model 



In the Glauber model [13], the starting point is the Woods-Saxon density distribution for 
nuclei, 

Po 



Pa(x) 



1 + exp [(|x| - Ro)/x] 



(4.26) 



where for a gold nucleus with weight A = 197 we use i?o = 6.4 fm and x = 0-^4 fm. The 
parameter po is chosen such that J (i^x/9^(x) = A. One can then define the nuclear thickness 
function 



dzpA{^) 



(4.27) 



and subsequently the number density of nucleons participating in the collision (npart) and 
the number density of binary collisions (ncoii)- For a collision of two nuclei with weight A 



^Note that a small numerical error occurred in the linearized hydrodynamic solver and the corresponding 
figure in Ref. [4^. This error has been corrected in [Figure 4.4[ 
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at an impact parameter b, one has 



npi,rt{x, y,b) = Ta \ X + -,y 



+ Ta{ X - -,y 



f^Coii(a;, y, h) = aTA ( x + ^, y j Ta f x - ^, y 



(tTa ix 




(4.28) 



where a is the nucleon-nucleon cross section. We assume a ~ 40 mb for Au+Au cohisions 
at y/s = 200 GeV per nucleon pair. 

While the total number of participating nucleons A^part(^) = / dxdynpi^j.t{x , y, b) will be 
used to characterize the centrality class of the collision, as an initial condition for the energy 
density we will only use the parametrization 



e(r = To, x, y, b) = const x ncoll(2;, y, i 



(4.29) 



since it gives a sensible description of the multiplicity distribution of experimental data, 
as will be discussed later on. In the following, "Glauber-model initial condition" is used 



synonymous to Equation 4.29 



The constant in Equation 4.29| is chosen such that the central energy density for zero 
impact parameter, e(r = tq, 0,0,0) corresponds to a predefined temperature Tj via the 
equation of state. This temperature will be treated as a free parameter and is eventually 
fixed by matching to experimental data on the multiplicity. 



4.3.2 The CGC model 

The other model commonly used to obtain initial conditions for hydrodynamics is the so- 
called Color-Glass-Condensate approach, based on ideas of gluon saturation at high energies. 
In particular, we use a modified version of the KLN (Kharzeev-Levin-Nardi) /cy-factorization 
approach [sil, due to Drescher et al. [S]. We follow exactly the procedure described in 71 1 
and in fact we use the same numerical code, provided to us by the authors and only slightly 
modified to output initial conditions suitable for input into our viscous hydrodynamics 
program. 
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(P^TdY J pI 



In this model, the number density of gluons produced in a colhsion of two nuclei with 
atomic weight A is given by 

— — / d^kT as(A:T) (pT + kT)^/4;xT) (/'a(2;2, (pT-kT)^/4;xT) 

(4.30) 

where pt and Y are the transverse momentum and rapidity of the produced gluons, respec- 
tively. xi^2 = Pt ^ exp(iby)/y^ is the momentum fraction of the colliding gluon ladders 
with ^/s the center of mass collision energy and as{kT) is the strong coupling constant at 
momentum scale = \^t\- 

The value of the normalization constant M is unimportant here, since as for Glauber 
initial conditions, we treat the overall normalization of the initial energy density distribution 
as a free parameter. The unintegrated gluon distribution functions are taken as 



1 Ql 

UsiQl) max((5s, krp 



(x,4;xt) = — ^ ;^^-2-P(xt)(1-x)4 , (4.31) 



P{^t) is the probability of finding at least one nucleon at transverse position x-p, taken 
from the definition for npart 

nxT) = l-(l-^)^ (4.32) 

where Ta and a are as defined in the previous section. 

The saturation scale at a given momentum fraction x and transverse coordinate x^ is 
given by 

,,.,..,..CeV^(5i™)(M)\ 

The growth speed is taken to be A = 0.288. 

The initial conditions for hydrodynamic evolution require that we specify the energy 
density in the transverse plane at some initial proper time tq at which the medium has 
thermalized. Equation 4.30[ on the other hand, is in principle valid at a time Tg = 1/Qs 
at which the medium is likely not yet in thermal equilibrium. To obtain the desired initial 



conditions, we again follow 



711] and assume that the number of gluons is effectively con- 



served during the evolution from to tq and so the number density profile is the same at 
both times, scaled by the one-dimensional Bjorken expansion n(ro) = ^n(rs). The energy 
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density can then be obtained from the number density through thermodynamic relations — 
it is proportional to the number density to the 4/3 power. Again, we take the overall 
normalization as a free parameter, so the initial energy density is finally given as 



e(r = To, XT, h) = const x 



:(XT,6) 



n4/3 



(4.34) 



jPxTdY 

where the number density is given by Equation 4.30| evaluated at central rapidity 1" = 0. 



As a final comment, it shou 
McLerran-Venugopalan model 
be discussed in the next-section. 



d bepointed out that the original version of the CGC, the 
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94( 1 ■ differs from the KLN ansatz we used here, as will 



4-3.3 Spatial and momentum anisotropy 

One of the key parameters discussed in the following is the eccentricity (or spatial anisotropy) 
of the collision geometry. Following [3] , we define it as 



(4.35) 



where {)e denotes an averaging procedure over space with the energy density e as a weighting 
factor. Shown in Figure CS^ , a plot of for different centralities highlights the quanti- 
tative difference between the initial energy density from the Glauber and CGC model. 
Equation 4.29| and Equation 4.34[ respectively. As can be seen from this figure, the CGC 
model generally gives a higher spatial anisotropy than the Glauber model. Note that 
the results for the CGC model shown here are extreme in the sense that the McLerran- 
Venugopalan model gives spatial eccentricities which essentially match the ones from the 
Glauber model [95[. This allows us to use the difference between the CGC and Glauber 



models as an indication of the systematic theoretical error stemming from the ignorance of 
the correct physical initial condition. 

Hydrodynamics converts pressure gradients into fluid velocities, and hence one expects 
the spatial anisotropy to decrease at the expense of a momentum anisotropy (which is related 
to the magnitude of the elliptic flow). We follow 96| in defining a momentum anisotropy 
according to 



(rji'xx _j_ J-j/y) ' 



(4.36) 
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(a) (b) 




Figure 4.5: Left: The initial spatial anisotropy for the Glauber and CGC model. Right: 
The time evolution of the spatial and momentum anisotropy for a collision with 5 = 7 fm 
in ideal hydrodynamics. 



where we stress that here () denotes spatial averaging with weight factor unity. Figure 4T5| 3 



shows the time evolution in ideal hydrodynamics (r//s <^ 1) of both the spatial and momen- 
tum anisotropies for a heavy- ion collision at 6 = 7 fm modeled through Glauber and CGC 
initial conditions. As one can see, for the same impact parameter, the higher initial spatial 
anisotropy for the CGC model eventually leads to a higher momentum anisotropy than the 
Glauber model. Using a quasiparticle interpretation where the energy momentum tensor is 
given by 

the momentum anisotropy Cp can be approximately related to the integrated elliptic flow 



96|, 



97|. We find this proportionality to be 



v^2^(j3), with a proportionality factor of ~ 2 
maintained even for non- vanishing shear viscosity, as can be seen in [Figure 49} 

4.4 Results 

4-4- 1 Which parameters matter? 

In the following, we will attempt to obtain limits on the mean value (throughout the hy- 
drodynamic evolution) of the ratio rj/s from experimental data. While e.g. temperature 
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Figure 4.6: Spatial and momentum anisotropy for the Glauber model at 6 = 7 fm with 
Tj = 0.353 GeV, tq = 1 fm/c and various values for the viscosity (grid spacing a = 2 GeV~^). 
(a): The dependence on the initialization value of the shear tensor: shown are results for 
vanishing initial value (n|^'^^ = 0) and Navier-Stokes initial value (J^i^^^ 7^ 0), given in 
Equation 4.4[ (b): The dependence on the choice of value for Tt^,Xi: shown are results 
for = Ai = (labelled "IS") and = ^i2^^R^ ^1 = 2^ (labelled "AdS"). For 
= ^j^^ results for Ai = (not shown) would be indistinguishable by bare eye 

from those for Ai = 
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variations of r]/ s are to be expected in the real physical systems, probing for such variations 
would invariable force us to introduce more unknown parameters. We prefer to leave this 
program for future studies once robust results for the mean value of r]/s exist. Having fixed 
the equation of state and the freeze-out procedure as explained in the previous sections, the 
remaining choices that have to be made in the hydrodynamic model are the 

• Initial energy density profile: Glauber or CGC 

• Initial value of shear tensor: vanishing or Navier-Stokes value 

• Hydrodynamic starting time tq 

• Second-order coefficients: relaxation time r^r and Ai 



Ansatz for non-equilibrium particle distribution Equation 4.14 



where it is to be understood that we fix the initial energy density normalization (Tj) and 
the freeze-out temperature Tf such that the model provides a reasonable description of 
the experimental data on multiplicity and {pt)- Historically, a strong emphasis has been 



99l |. For this reason, we 



put on requiring a small value of tq for ideal hydrodynamics 
will discuss the dependence on tq separately in [section 4.4.4I A good indicator for which 
parameters matter is the momentum anisotropy since it is very sensitive to the value of 



ri/s. From Figure 4.5 one therefore immediately concludes that the choice of Glauber or 
CGC initial conditions is important since it has a large effect on Cp. Fortunately, most of 
the other choices turn out not to have a strong influence on the resulting V2 coefficient and 
hence the extracted rj/s. In the following we test for this sensitivity by studying Cp for a 
"generic" heavy-ion collision of two gold nuclei, modeled by Glauber initial conditions at 
an initial starting temperature of Tj = 0.353, an impact parameter of 6 = 7 fm, and various 
choices of the above parameters. 



Figure 4.6 shows the time evolution of e^, Cp for various values of r]/s. From these plots, 
it can be seen that Cp (and hence V2) clearly is sensitive to the value of rj/s, suggesting 
that it can be a useful observable to determine the viscosity of the fluid from experiment. 
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Figure 4.7: Charged hadron elliptic flow for the Glauber model at 6 = 7 fm with Tj = 0.353 
GeV, To = 1 fm/c and various viscosities. 



However, in order to be a useful probe of the fluid viscosity, the dependence of the final 
value of Bp on other parameters should be much weaker than the dependence on rj/s. In 



Figure 4.6 1 we show Cp, calculated for U^'^^tq) = and 11^'^ (tq) equal to the Navier- 



Stokes value, lEquation 4.41 As can be seen from this figure, the resulting anisotropics are 


essentially independent of this choice, corroborating the finding in Ref. [551. 157l|. Similarly, 



m 



Figure CBj j we show Cp calculated in simulations where the values of the second-order 



transport coefficients were either those of a weakly-coupled Miiller-Israel-Stewart theory 
(ttt = 6^, Ai = 0) or those inspired by a strongly coupled M = 4 SYM plasma (r^ = 
2(2 — In 2)^, Ai = 2^)- Again, the dependence of ep on the choice of the values of 
r,r,Ai can be seen to be very weak for the values of i]/s shown here. This result is in 
stark contrast to the findings of Ref. [55[], where a large sensitivity on the value of was 



found. However, recall that Ref. 



55l | used evolution equations that differ from Equation 4.1 



and in particular do not respect conformal invariance. As argued in [section 4.2.21 it is 
therefore expected to encounter anomalously large sensitivity to the value of the second 
order transport coefficients. 

To study the dependence of results on the ansatz of the non-equilibrium particle distribu- 
tion function ()4.14p , one would want to quantify the effect of neglecting terms of higher order 
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in momenta in Equation 4.14 To estimate this, let us rewrite Ed? N / d^p = Ed'^N^^^ /d^p + 



Ed'^N^^^ /d^p, where A^^'^) contains only the equilibrium part where f{x^,p^) = Jq ^ 



and perform a Pade-type resummation, 

^^SjyPade ^ ^^3^(0) ^ 

d^p d^p I 



- ^ -.3- . .3^(.) ,3p • (4.38) 



Since Equation 4.38 contains powers of momenta to all orders when re-expanded, the differ- 



ence between the ansatz (|4.14p and the Pade resummed particle spectra can give a handle 



on the systematic error of the truncation used in Equation 4.14 Shown in Figure 4.7, this 
difference suggests that this systematic error is small for momenta pT ^ 2.5 GeV. There- 
fore, we do not expect our results to have a large systematic uncertainty coming from the 
particular ansatz (j4.14p for these momenta. 

To summarize, for values of r]/s < 0.2, the results for the momentum anisotropy are 
essentially insensitive to the choices for the second-order transport coefficients Tt^,\i and 
the initialization of the shear tensor 11'^'^ (r = tq). Conversely, Cp is sensitive to the value 
of viscosity and the choice of initial energy density profile (initial eccentricity). Since the 
physical initial condition is currently unknown, this dependence will turn out to be the 
dominant systematic uncertainty in determining rj/ s from experimental data. 

4- 4-2 Multiplicity and radial flow 

As outlined in the introduction, we want to match the hydrodynamic model to experimental 
data for the multiplicity, thereby fixing the constant in Equations 14.291 and 14.341 This 
translates to fixing an initial central temperature Tj for 6 = 0, which we will quote in the 
following. 

For a constant speed of sound, the evolution for ideal hydrodynamics is isentropic, 
while for viscous hydrodynamics additional entropy is produced. Since the multiplicity is 
a measure of the entropy of the system, one expects an increase of multiplicity for viscous 
compared to ideal hydrodynamic evolution. This increase in final multiplicity has been 
measured as a function of rj/ s for the semi-realistic speed of sound [Figure -O] in central 



heavy-ion collisions in Ref. 



47[, and found to be approximateljo a factor of 0.75r//s. (See 



The quoted fraction is for a hydrodynamic starting time of ro = 1 fm/c. Reducing ro leads to consid- 
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Ref. 



71 



loot] for related calculations in simplified models.) Reducing Ti accordingly there- 
fore ensures that for viscous hydrodynamics, the multiphcity in central colhsions will stay 
close to that of ideal hydrodynamics. 

Hydrodynamics gradually converts pressure gradients into flow velocities, which in turn 
relate to the mean particle momenta. Starting at a predefined time tq and requiring the 
hydrodynamic model spectra to match the experimental data on particle {pt) then fixes the 
freeze-out temperature Tf. 

For both Glauber-type and CGC-type model initial conditions, the experimental impact 
parameter dependence of the multiplicity and {px) is reasonably well parametrized for both 
ideal hydrodynamics as well as viscous hydrodynamics provided Tj is adjusted accordingly 
(see [Figure 48] ) . The values for T, used in the simulations are compiled in ITable 4!T1 We 
recall that no chemical potential is included in our equation of state, prohibiting a distinction 
between particles and anti-particles, and chemical and kinetic freeze-out of particles occurs 
at the same temperature. Furthermore, approximating the equilibrium particle-distributions 
for bosons by a Boltzmann distribution (j4.14p leads to small, but consistent underestimation 
of the multiplicity of light particles, such as pions. For these reasons, it does not make sense 
to attempt a precision fit to experimental data, especially for pions and protons. Rather, 
we have aimed for a sensible description of the overall centrality dependence of multiplicity 
and {pt) of kaons. 

Note that in particular for the CGC model one could achieve a better fit to the data 
on mean {pt) by increasing the freeze-out temperature by ~ 10 MeV. This would also lead 
to a decrease in elliptic flow for this model. However, to facilitate comparison between the 
CGC and Glauber initial conditions, we have kept Tf the same for both models. 



^.4-3 Elliptic flow 

Having fixed the parameters tq, Ti,Tf for a given r//s to provide a reasonable description of 
the experimental data, a sensible comparison between the model and experimental results 
for the elliptic flow coefficient can be attempted. For charged hadrons, the integrated and 



erably larger entropy production. 
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Figure 4.8: Centrality depen dence of total multiplicity dN/dY and {pt) for vr"'", vr , , 



K-, p and p from PHENIX 



lOll ] for Au+Au collisions at ^/s = 200 GeV, compared to 



the viscous hydrodynamic model and various rj/s, for Glauber initial conditions and CGC 
initial conditions. The model parameters used here are tq = 1 fm/c, = Qr]/s, Ai = 0, 
Tf = 140 MeV and adjusted Ti fsee lTable 4.T1) . 
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Table 4.1: Summary of parameters used for the viscous hydrodynamics simulations 





'11 * 


± J [Kjl£ V J 


If [»jevj 




a [GeV-^l 
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Glauber 


10-4 
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minimum-bias V2 coefficients are shown in Figure 4.9 for Glauber and CGC initial conditions. 
As noted in [section 4.3.31 charged hadron u™* turns out to be very well reproduced by the 
momentum eccentricity ^ Cp, evaluated when the last fluid cell has cooled below Tf. This 
agreement is independent from impact parameter or viscosity and hence may serve as a 
more direct method on obtaining an estimate for w™' if one cannot (or does not want to) 
make use of the Cooper- Frye freeze-out procedure described in [section 4.2.7[ 

The comparison of the hydrodynamic rnodel to experimental data with 90% confidence 
level systematic error bars from PHOBOS 102l | for the integrated elliptic flow in Figure 4.9 



suggests a maximum value of r//s ~ 0.16 for Glauber-type and r]/s ~ 0.24 for CGC-type 
initial conditions. Whereas for Glauber initial conditions, ideal hydrodynamics {r]/s ~ 0) 
gives results consistent with PHOBOS data, for CGC initial conditions zero viscosity does 



not give a good fit to the data, which is consistent with previous findings [85]. 

For minimum-bias f2, to date only experimental data using the event-plane method 
are available, where the statistical, but not the systematic error of that measurement is 
directly accessible. The dom inan t source of systematic error is associated with the presence 
of so-called non-flow effects 1051 ] . Recent results from STAR suggest that removal of these 



non -flow effects imply a reduction of the event-plane minimum bias V2 by 20 percent 



103 



lOJ]. For charged hadrons, a comparison of both the event-plane and the estimated non- 
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Figure 4.9: Comparison of hydrodynamic models to experimental d ata o n charged hadron 



integrated (left) and minimum bias (right) elliptic flow by PHOBOS 102t | and STAR 103l | 



respectively. STAR event plane dat a has been reduced by 20 percent to estimate the removal 
of non-flow contributions 



1031, 



lOJ]. The line thickness for the hydrodynamic model curves 
is an estimate of the accumulated numerical error (due to, e.g., finite grid spacing). The 
integrated V2 coefficient from the hydrodynamic models (full lines) is well reproduced by 
(dots); indeed, the difference between the full lines and dots gives an estimate of the 
systematic uncertainty of the freeze-out prescription. 
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flow corrected experimental data from STAR with the hydrodynamic model is shown in 
Figure 4.9 

For Glauber-type initial conditions, the data on minimum-bias V2 for charged hadrons 
is consistent with the hydrodynamic model for viscosities in the range rj/s G [0,0.1], while 
for the CGC case the respective range is rj/s £ [0.08,0.2]. It is interesting to note that for 
Glauber-type initial conditions, experimental data for both the integrated as well as the 
minimum-bias elliptic flow coefficient (corrected for non-flow efl'ects) seem to be reproduced 
besJzlby a hydrodynamic model with rj/s = 0.08 ~ This number has flrst appeared in the 
gauge/string duality context [3Q] and has been conjectured to be the universal lower bound 
on r]/s for any quantum field theory at finite temperature and zero chemical potential \w^. 
For CGC-type initial conditions, the charged hadron V2 data seems to favor a hydrodynamic 
model with r]/s ~ 0.16, well above this bound. 



4-4-4 Early vs. late thermalization 

Currently, there seems to be a common misunderstanding in the heavy-ion community that 
hydrodynamic models can universally only reproduce experimental data if they are initial- 
ized at early times tq < 1 fm/c. This notion has been labeled "early thermalization" and 
continues to create a lot of confusion. In this section, we argue that the matching of hydro- 
dynamics to data itself does not require tq < 1 fm/c. It is the additional assumptions about 
pre-equilibrium dynamics that lead to this conclusion for the Glauber initial conditions. 

Performing hydrodynamic simulations in the way we have described earlier, the energy 
density distribution is specified by either the Glauber or CGC model at an initial time tq. In 



Figure 4.10 we show the result for the elliptic flow coefficient (or the momentum anisotropy) 
for three different values of tq, namely 0.5, 1 and 2 fm/c, where also Tj and Tf have been 
changed in order to obtain roughly the same multiplicity and mean pT for each tq. As can 
be seen from this figure, the resulting final elliptic flow coefficient is essentially independent 
of the choice of tq. In particular, this implies that experimental data for bulk quantities 

^ In Ref. [isj a lower value of rj/s for the Glauber model was reported. The results for viscous hydrody- 
namics shown in [Figure 4.9| are identical to Ref. [i^, but the new STAR data with non-flow corrections 
became available only after had been published. 
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(a) (b) 




Figure 4.10: Momentum anisotropy (a) and elliptic flow for charged hadrons (b) for 6 = 7 
fm, r]/s = 0.08 and different hydrodynamic initialization times tq. Horizontal light gray 
lines in (a) are visual aids to compare the final value of Cp. As can be seen from these 
plots, neither the final Cp nor the charged hadron V2 depend sensitively on the value of tq 
if the same energy distribution is used as initial condition at the respective initialization 
times. Simulation parameters were Tj = 0.29 GeV, Tf = 0.14 GeV for tq = 2 fm/c, 
Ti = 0.36 GeV, Tf = 0.15 GeV for tq = 1 fm/c, and Tj = 0.43 GeV, Tf = 0.16 GeV for 
To = 0.5 fm/c. 
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can be reproduced by hydrodynamic models also for large initialization times, so no early 
thermalization assumption is needed. 

However, it is true that the above procedure assumes that the energy density distribu- 
tion remains unchanged up to the starting time of hydrodynamics, which arguably becomes 



93] to mimic the 



increasingly inaccurate for larger tq. It has therefore been suggested 
pre-hydro time evolution of the energy density distribution by assuming free-streaming 
of partons. Assuming free-streaming gives the maximal contrast to assuming hydrody- 
namic evolution, since the latter corresponds to very strong interactions while the former 
corresponds to no parton interactions at all. Indeed, one can calculate the effect of the 
free-streaming evolution on the spatial anisotropy, finding 98 1 



e.(r) = ^^, {R')= ^f^r"^ • (4-39) 

This implies that the spatial anisotropy decreases with time, whereas one can show that 
free-streaming does not lead to a build-up of e^. In other words, the eccentricity gets 
diluted without producing elliptic flow, such that once hydrodynamic evolution starts, it 
will not lead to as much V2 as it would have without the dilution effedo. It is tempting 
to conclude from this that by comparing to experimental data on elliptic flow one could 
place an upper bound on the maximally allowed dilution time, and interpret this as the 
thermalization time of the system. One should be aware, however, that this bound will 
depend on the assumption made about the pre-hydro evolution. Furthermore, one should 
take into account the fact that the initial state of the system remains unknown. For instance, 
the system could start with an energy density distribution similar to the CGC model, which 



has a fairly large eccentricity. Figure 4.11| shows that when allowing the eccentricity to get 



diluted according to Equation 4.39, it takes a time of r ~ 1.5 fm/c until the eccentricity 
has shrunk to that of the Glauber model. This implies that even when assuming no particle 
interactions (no elliptic flow build-up) for the first stage of the system evolution, one can 

* It seems that if one forces the energy-momentum tensor at the end of free-streaming period to match to 
that of ideal hydrodynamics (instantaneous thermahzation ), the res ulting fluid velocities are anisotropic, 
i.e. correspond to a non- vanishing elliptic flow coefficient [IOTI . llOSi ]. It is possible that this effect stems 
from neglecting velocity gradients (viscous hydrodynamic corrections) in the matching process. We ignore 
the complications of the detailed matching from free-streaming to hydrodynamics in the following. 
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Figure 4.11: Spatial eccentricity for the Glauber and CGC model compared to evolving the 
CGC model according to Equation 4.39| for r = 1.5 fm/c. This implies that starting with 
Glauber-type initial conditions at tq > 1 fm/c may not be unreasonable. 



get eccentricities which are Glauber-like after waiting for a significant fraction of the system 
life time. Allowing at least some particle interactions (which is probably more realistic), one 
expects some build-up of elliptic flow already in the dilution (or pre-equilibrium) phase, and 
therefore dilution (or "thermalization" ) times of r ~ 2 fm/c seem not to be incompatible 
with the observed final elliptic fiow even for non- vanishing viscosity. 

4.5 Summary and Conclusions 

We applied conformal relativistic viscous hydrodynamics to simulate Au+Au collisions at 
RHIC at energies of y/s = 200 GeV per nucleon pair. Besides one first-order transport 
coefficient (the shear viscosity) in general there are five second-order transport coefficients 
in this theory, for which one would have to supply values. We provided arguments that 
physical observables in the parameter range accessible to hydrodynamics (low momenta, 
central to semi-central collisions) do not seem to be strongly dependent on specific (rea- 
sonable) choices for any of these second-order coefficients. On the other hand, we do find 
a pronounced dependence of the elliptic fiow coefficient on the ratio of shear viscosity over 
entropy density, which suggests that by combining viscous hydrodynamics and experimen- 
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tal data a measurement of the quark-gluon plasma viscosity may not be futile. However, 
we have shown that our ignorance about the precise distribution of energy density at the 
earliest stages of a heavy-ion collision introduces a large systematic uncertainty in the final 
elliptic flow of the hydrodynamic model. Adding to this is the considerable experimental 
uncertainty pertaining to the removing of non-flow contributions to the elliptic flow. For 
these reasons, we are unable to make precise statements about the value of the shear vis- 
cosity of the quark-gluon plasma and in particular cannot place a firm lower bound on r//s. 
Indeed, our hydrodynamic models seem to be able to consistently describe experimental 
data for multiplicity, radial flow and elliptic flow of bulk charged hadrons for a wide range 
of viscosity over entropy ratios. 



7] 



0.1 lb O.l(theory) it 0.08(experiment), 



(4.40) 



where we estimated the systematic uncertainties for both theory and experiment from the 
results shown in Figure 4.9 , We stress that Equation 4.40| does not account for physics not 
included in our model, such as finite chemical potential, bulk viscosity, heat flow, hadron 
cascades, three-dimensional fluid dynamic effects and possibly many more. Consistent in- 
clusion of all these may result in changes of the central value and theory uncertainty in 
[Equation 4.40[ Nevertheless, none of the mentioned refinements is currently expected to 
dramatically increase the elliptic flow coefficient (th ough some increase may be expected 



when e.g. implementing partial chemical equilibrium 



1091]). Therefore, we seem to be able 



to exclude viscosities oi rj/s > 0.5 with high confldence, which indicates that the quark- 



gluon plasma displays less friction than a ny other know n laboratory fluid 



groups have come to similar conclusions 
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106 



lld |. Other 



113|. 



To better quantify the shear viscosity of the quark-gluon plasma at RHIC calls for 
more work, both in theory and experiment. On the theory side, a promising ro ute s eems 



to 



3 6 the study of fluct u ation s and comparing to existing experimental data jl02 
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ILl 



llfil . 



Ill Ill8l . 



119|. For instance, it might be interesting to investigate the 



critical value oi rj/s for the onset of turbulen ce in heavy- ion collisions and explore possible 
consequences of fully developed turbulence |;120]. However, maybe most importantly, a 
more thorough understanding of the earliest stages of a heavy-ion collision, in particular 
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thermalization, could fix the initial conditions for hydrodynamics and hence dramatically 
reduce the theoretical uncertainty in final observables. 

Leaving these ideas for future work, we stress that with the advent of conformal rel- 
ativistic viscous hydrodynamics at least the uncertainties of the hydrodynamic evolution 
itself now seem to be under control. We hope that this serves as another step towards a 
better understanding of the dynamics of relativistic heavy-ion collisions. 
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Chapter 5 

VISCOUS HYDRODYNAMIC PREDICTIONS FOR THE LHC 
5.1 Introduction 



Using the knowledge gained from viscous hydrodynamic simulations for RHIC, it should 
be possible to predict experimental results at the Large Hadron Collider (LHC), which will 
collide lead ions at a maximum center of mass energy of ^/s = 5.5 TeV per nucleon pair 
compared to ^/s = 200 GeV gold ions at RHIC. If experimental data on, e.g., V2 from LHC 
is close to the hydrodynamic model prediction, this would confirm that real progress has 
been made in understanding nuclear matter at extreme energy densities; if far away, it may 
indicate that the successful hydrodynamic description of experimental data from RHIC was 
a coincidence. 

Regardless of the outcome, the advent of the RHIC experiments clearly has lead to 
major progress in the theory and application of hydrodynamics to heavy- ion collisions. A 
few years ago the form of the hydrodynamic equations in the presence of shear viscosity 
71 was s till unres olved, with different groups keeping some terms while neglecting others 
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121 



122l |. For the case of approximately conformal theories, where the viscosity 



coefficient for bulk — but not shear — becomes negligible, all possible terms to second order in 



gradients were derived in Ref. [3l|] (see chapter 3 ) , and their relative importance investigated 
in Ref. [l| ( [chapter"! ) . Thre e of the g roups performing viscous hydrodynamic simulations 



now agree on these terms 



[],[i3, 



124l |. while another group 



uses a different formalism 



that gives matching results. While this development still leaves out the consistent treatment 
of bulk viscosity, the quantitative suppression of elliptic flow by shear viscosity is therefore 
essentially understood. As show n in 



simulations to experimental data 



lOl 



chapter 4 from comparison of viscous hydrodynamic 



103( 1 ■ one can infer an upper limit of the ratio of shear 



viscosity over entropy density, ij/s < 0.5, for the matter produced in Au+ Au cqllisions a t 



s = 200 GeV [li, which is in agreement to extractions by other methods 



111 



112 



113] 
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A sizeable uncertainty for this limit comes from the fact that the initial conditions for 
the hydrodynamic evokition are poorly known, with the two main models, the Glauber 
and Color-Glass-Condensate (CGC) models, giving different results for the elliptic flow 
coefficient This difference can be understood to originate from the different initial 
spatial eccentricity tx in the Glauber/CGC models, which we recall is defined as 



(y2 - X 



2\ 

^, (5.1) 



(y2 + x'^)^ ' 

where the symbols ()^ denote averaging over the initial energy density in the transverse 
plane, e(x, y). 

Indeed, it had been suggested 125l | that the elliptic flow coefficient V2 at the end of the 
hydrodynamic evolution would be strictly proportional to the initial spatial eccentricity, 
v^l^x oc const., if the fluid was evolving without any viscous stresses for an infinitely long 
time. This is to be contrasted with experimental data indicat ing a proportionality factor 



of total multiplicity over overlap area V2lex oc dN/ dY/ ^overlap 125l | . Total multiplicity '^^ 



dY 



here refers to the total number of observed particles N per unit rapidity Y, while the overlap 
area is calculated as 



•^overlap = T^V {x'^) iv'^) ■ (5.2) 

Ideal fluid dynamics does not adequately describe the latest stage of a heavy-ion collision 
(the hadron gas), because of the large viscosity coefficient in this stage |l2^. Therefore, 
the hydrodynamic stage lasts only for a finite time (e.g., until all fluid cells have cooled 
below the decoupling temperature), resulting in a dependence of ^2/62; on dN/dY. Also, 
viscous effects affect the proportionality between V2 and Cx, leading to a behavior that is 
qualitatively similar to that observed in the data [123.]. 

One of the objectives of this work is to extend the energy range for fluid dynamic results 
of V2/ex from Au-|-Au collisions at top RHIC to Pb-|-Pb collisions at top LHC energies, 
as well as to study the dependence on shear viscosity. If in the future either or the 
mean ij/s becomes known, these results can thus be used to constrain the respective other 
quantity from experimental data. On the other hand, the values of shear viscosity for 
which the Glauber/CGC models match to experimental data at top RHIC energies have 
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been extracted in Ref. for Au+Au collisions. Since rj/s averaged over the system 

evolution is not expected to be dramatically different for Pb+Pb collisions at the LHC, 
another objective of this work is to obtain a prediction for the elliptic flow coefficient for 
the LHC based on the best-fit values to RHIC. 

Finally, the feasibility of det ecting elliptic flow in p+p collisions at ^/s = 14 TeV at 



the LHC is being discussed jl27l ]. As a reference for other approaches and experiment, it 
interesting to study the possible size and viscosity dependence of V2 under the hypothetical 
assumption that the bulk evolution following p+p collisions could be captured by fluid 
dynamics. 

5.2 Setup 

To make predictions for nuclear collisions at LHC energies,^we use our hydrodynamic model 
that successfully described experimental data at RHIC [ll,!^ and make modifications to the 
input parameters appropriate for the higher collision energies at the LHC. 

As a reminder, the hydrodynamic model is based on the conservation of the energy 



momentum tensor 



4 



A^A^DH"^ + -n'^'^(V„n") 



Irj-^ 27] 2 

where e,p and are the energy density, pressure, and fluid 4- velocity, respectively. D = 
u^D^ and Va = AqZ)^ are time-like and space-like projections of the covariant derivative 
D^j, where A^'^ = g^'^ — u^u'^ and we remind the compact notations 

Ai^^B^) = (a^A(^ + A^A(^ - |A"/5a^^) and u^^ = ^ (y„n,,, - V^n^). For relativis- 



tic nuclear collisions it is convenient to follow Bjorken [128l | and use Milne coordinates 
proper time r = \/t'^ — and spacetime rapidity ^ = atanh|, in which the metric becomes 
9iLv = diag(l, —1, —1, — r^), and assume that close to ^ = 0, the hydrodynamic degrees of 
freedom are approximately boost-invariant ~ y). 

The hydrodynamic equations Dfj^T^'^ = then constitute an initial value problem in 
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Table 5.1: Central collision parameters used for the viscous hydrodynamics simulations 
{Tf = 0.14 GeV for all). 



Beam 


Initial cond. 


dY 


Ti [GeV] 


^ [GeV] 


To [fm/c] 


Gold 


Glauber 


800 


0.34 


200 


1 


Gold 


CGC 


800 


0.31 


200 


1 


Lead 


Glauber 


1800 


0.42 


5500 


1 


Lead 


CGC 


1800 


0.39 


5500 


1 


Protons 


Glauber 


6 


0.400 


14000 


0.5 


Protons 


Glauber 


6 


0.305 


14000 


1 


Protons 


Glauber 


6 


0.270 


14000 


2 



n 

proper time and transverse space, and are solved numerically (see Ref. [1]). The input 
parameters for hydrodynamic evolution are the equation of state p = p(e) and the first 
(second) order hydrodynamic transport coefficients rj (r,ri Ai, A2, A3). The values for Ai^2,3 
have been found to hardly affect the boost-invariant hydrodynamic evolution for Au+Au 
collisions at RHIC [1], so here they are generally set to zero. 

The mechanisms leading to thermalization (the onset of hydrodynamic behavior) are 
not well understood in nuclear collisions. Therefore, it is not known how the thermalization 
time To at which hydrodynamic evolution is started depends on the collision energy. Barring 
further insight, we start hydrodynamic evolution for the LHC at the same time as for the 
RHIC simulations (tq = 1 fm/c). At this time, the initial conditions for the transverse energy 
density e{x,y) are given by the Glauber or CGC model, respectively, the fluid velocities are 
assumed to vanish, and the shear tensor 11^'' is set to zero (other values for 11'^'^ do not 
seem to affect the final results [l|, [s^] ) • For brevity, we refer to [chapter ~4: for the details 
of the Glauber and CGC models, but note that we use the Woods-Saxon parameters of 
radius Rq = 6.4(6.6) fm and skin depth x = 0.54(0.55) fm for gold (lead), and assume a 
nucleon-nucleon cross section of o" = 40 (60) mb for ^/s = 200 (5500) GeV collisions. 

The overall normalization of the initial energy density (parametrized by the initial tem- 



64 



perature at the center Tj) was adjusted to match the experimentally observed multiplicity at 



rh: 

ity 



C; by analogy, for LHC the normalization is adjusted to match the predicted multiplic- 



13C 
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132]. Since we lack detailed knowledge about its temperature dependence, 



the ratio of shear viscosity to entropy density ij/s is set to be constant during the hydro- 
dynamic evolution (equal to the average ove r th e spacetime evolution of the system). The 



relaxation time coefficient r^r is expected 
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1331 ] to lie in the range ^(e+p) — 2.6 — 6. The 



equation of state (EoS) can in principle be provided by lattice QCD. While at present there 
are points of disagreement between lattice groups about, e.g., the preci se locatio n of the 
QCD phase transition, there is consensus that it is an analytic crossover jl34l . Il35| ]. There- 
fore, we use a lattice-inspired EoS [3] that is consistent with both the current consensus 



1351 ] ■ we expect that using a different lattice 



and perturbative QCD; also, since it resembles 
EoS will have a minor effect on our results. In fact, as a preview of work in progress, some of 
the calculations were also done with an equation of state taken from Ref. 



1351] ( [Figure 5.5 



at the end of the chapter) and the results compared (see Figure 5.2). 

Once a given fluid cell has cooled down to the decoupling temperature Tj, its energy and 
momentum are converted into particle degrees of freedom using the Cooper-Frye freeze-out 
prescription [t^. A value of Tj = 0.14 GeV was determined by matching to RHIC data and 
will also be used for LHC energies, assuming that it is mostly determined by local conditions, 
and less so by initial energy density, system size or collision energy. The distribution of the 
particle ^ degrees of freedom may be further evolved using a hadronic cascade code (as in 
Ref. [83i]), or in a more simple approach the unstable particle resonances are allowed to 
decay, without further evolving the stable particle distributions. In both cases, the total 
multiplicity and particle correlations (such as the elliptic flow coefficient) are then calculated 
from the stable particle distribution (cf. [1]). Surprisingly, it was found in Ref. [l], 3] (see 



chapter 4) that the momentum integrated elliptic flow coefficient for charged hadrons — to 



good approximation — is equal to half the momentum anisotropy, 

1 I fdxdyT'''' -ryy 

2 ^ 2 / dxdy T^^ + Tw 



(5.3) 



Since the momentum anisotropy is a property of the fluid, it is independent on the details of 
the freeze-out procedure and only mildly dependent on the choices of tq, Ty. Unlike at RHIC 
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where pairs of tq and Tj could be fine-tuned to fit the particle spectra at central collisions, 



no such extra information is available for the LHC. Hence Equation 5.3 may provide the 
most reliable way of determining the elliptic flow of charged hadrons, and will be used in 
the following. Similarly, one can use the total entropy per unit spacetime rapidity in the 
fluid as a proxy for the total (cha rged hadron) multiplicity per unit rapidity (w^) with 



a proportionality factor 



13 



136 1 

dS dS 



~ 4.87 

dY dY 



7.85 



(5.4) 



d^ dY dY dY 

Note that for a gas of massive hadrons in thermal equilibrium at Tj = 0.14 GeV the ratio 
of entropy to particle density is ~ 6.41, but the decay of unstable resonances produces 



additional entropy, resulting in [Equation 5.4 Since results from RHIC suggest there is onl y 
approximately 10% viscous entropy production during the hydrodynamic phase 47, 123| . 



the entropy ^ at r = tq can be used to estimate the final particle multiplicity. In the 
case of the LHC, the world average f or th e predicted charged hadron multiplicity for central 



Pb+Pb collisions at ^/s = 5.5 TeV 



13Q|, 



dY 



1800, can be used to estimate the total 



entropy at r = tq , and hence the overall normalization Tj of the initial energy density (see 
ITable 5.ip . 

Using Equations 15.31 and 15.41 for the multiplicity and elliptic flow allows to make pre- 
dictions for the LHC without having to model the hadronic freeze-out, which should make 
the results more robust. However, as a consequence one does not get information about the 



momentum dependence of t he el 
predictions by other groups 
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iptic flow coefficient, prohibiting detailed comparison with 



137[. 



5.3 Results 



With the initial energy density distribution fixed at tq, the hydrodynamic model then gives 
predictions for the ratio of v^l&x at the LHC. In Figure h^\ the results are shown for 
three different values of shear viscosity, for two different initial conditions and two different 
beams/collision energies (Au+Au at \fs = 200 GeV, Pb-|-Pb at ^ = 5.5 TeV). The 
resulting values for seem to be quasi- universal functions of the total multiplicity 

scaled by the overlap area /Soveriap, only depending on the value oi rj/s (and, to a lesser 
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Figure 5.1: Anisotropy (j5.3p divided by (jS.ip . as a function of initial entropy ()5.4p divided 
by (15. 2p . Shown are results from hydrodynamic simulations for ^/s = 200 GeV Au+Au 
(RHIC) and ^/s = 5.5 TeV Pb+Pb collisions (LHC). For com parison, experimental data 



for V2 from RHIC 
measured 



dY 



101| 



138( 1 ■ divided by from two models |113l |. is shown as a function of 



divided by (|5.2p . See text for details. 



extent, the collision energy). The deviations of the RHIC simulations from the universal 
curve can be argued to arise from a combination of the finite lifetime of the hydrodynamic 
phase at ^/s = 200 GeV and the presence of the QCD ph ase transition, and is strongest for 
ideal hydrodynamics, in agreement with earlier findings 123l |. 



Also shown in Figure 5.1 is experimental data for the elliptic flow coefficient for Au+Au 
collisions at RHIC, normalized by Cx fr om a Monte-Carlo calculation (including fluctuations) 
in Glauber and CGC models (see Ref. |ll3l | for details). Since Cx is not directly measurable, 
the differently normalized data gives an estimate of the overall size of V2/ex at RHIC. 
Directly matching experimental data on V2 using a hydrodynamic model with an initial 
specified by the Glauber or CGC model, a reasonable fit was achieved for a mean value of 



67 



o 



0.08 - 

0.06 - 

0.04 - 

0.02 - 




Glauber, tj/s = 0.08 
CGC, Tj/s = 0.16 
PHOBOS (RHIC) 
I I I I 



_L 



_L 



50 100 150 200 250 300 350 400 
N 



Figure 5.2: Anisotropy (|5.3p prediction for ^/s = 5.5 TeV Pb+Pb collisions (LHC), as a 
function of centrality. Prediction is based on values of rj/s for the Glauber / CGC model 
that matched ^/s = 200 GeV Au+Au collision data from PHOBOS at RHIC ([l38|, shown 
for comparison). 



r]/s ~ 0.08 and rj/s ~ 0.16, respectively 

Under the assumption that the average rj/s is similar for collisions at RHIC and the 
LHC (along with the assumptions discussed above), one can make a prediction for the 
integrated elliptic flow coefficient for charged hadrons as a function of impact parameter 
(or more customarily the number of participants A'part; cf. l|). The result is shown in 
[Figure 5^ As can be seen, we expect integrated V2 at the LHC to be about ten percent 
arge r than at RHIC, which is less than the increase predicted b y ide al hydrodynamics 



139|, and in agreement with the extrapolations by Drescher et al. 13l|. For comparison. 



[Figure 5^ a) shows these LHC prediction curves along with those with rj/s set to 0.0001, 
corresponding to ideal hydrodynamics and illustrating the larger value of V2 predicted by 
ideal hydrodynamics. Also, as can be seen in Figure 5T3| b), remaining uncertainty in the 
equation of state seems to have little effect on this prediction. 
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Figure 5.3: Predicted anisotropy from [Figure 5.2[ in comparison to the value when 
T]/ s = 0.0001, corresponding to ideal hydrodynamics (a) (upper black curves), and using an 
alternate lattice QCD equation of state (b) (circles) — see [Figure 5.5 



Finally, using the charge density parametrization of the proton p{b) in Ref. [l4Cll | as an 
equivalent of the nuclear thickness function in the Glauber model (cf. [1]) one obtains an 
estimate for the shape of the transverse energy density foll owing a r elativistic p+p collision. 



Using the predicted multiplicity at mid-rapidity jy ~ 6 
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132] for ^ = 14 TeV p+p 



collisions at the LHC, one can again use Equation 5.4 to infer the overall normalization of 
the energy density (or Tj) at r = tq (see lTable 5.1|) . As a "Gedankenexperiment" one can 
then ask how much elliptic flow would be generated in LHC p+p collisions if the subsequent 
evolution was well approximated by boost-invariant viscous hydrodynamics. One finds that 
for ideal hydrodynamics ^ ~ f 2 ~ 0.035 for integrated \v2\ in minimum bias collisions 
(cf. (28) in l|), while for rj/s > 0.08, V2 typically changes by almost 100 percent when 



varying the relaxation time ^(e + p) between 2.6 and 6 and varying tq by a factor of two 



(see Figure 5.4). This indicates that for r]/s > 0.08, the hydrodynamic gradient expansion 
does not converge and as a consequence it is unlikely that elliptic flow develops in p+p 
collisions at top LHC energies. If experiments find a non-vanishing value for integrated 
1^2 1 > 0.02 in minimum bias p+p collisions, this would be an indication for an extremely 
small viscosity r]/s < 0.08 in deconfined nuclear matter. 
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Figure 5.4: Bands encompassing the calculated anisotropy curves for ^/s = 14 TeV p+p 
collisions. The relaxation time -^{e + p) was varied between 2.6 and 6 and tq from 0.5 fm 
to 2.0 fm for each value of rj/s. (Note that much of the 7]/s = 0.08 band is obscured by the 
Tj/s = 0.16 band, as both have significant dependence on the relaxation time and tq.) 



To conclude, viscous hydrodynamics can be used to make predictions for the ratio of 
V2I&X as a function of multiplicity and Assuming a multiplicity of ^^^^^ — 1800 for the 
matter created in Pb+Pb collisions at LHC, as well as ry/s similar to RHIC, we predict the 
integrated elliptic flow for charged hadrons to be 10% larger at the LHC than at RHIC. 
We expect V2 measurements in p+p collisions to be consistent with zero, unless the shear 
viscosity is extremely small (r//s < 0.08). 73| 
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Figure 5.5: Equation of state from Laine and Schr oder [3] — used for all the main results of 
chapters H] and [5] — and from Bazavov et al. ( \l3_^ , p4+HRG) — used to generate the circles 
Figure 5T3| ^b). Above are the energy density (a) and pressure (b) — scaled by — as a 



m 



function of temperature. Below (c) is the speed of sound squared. 
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Chapter 6 

FINAL STATE INTERACTIONS AND THE DISTORTED WAVE 

EMISSION FUNCTION 



Although much of the work in this chapter was performed before the viscous hydro- 
dynamic simulations described previously, it is most naturally described here, after some 
understanding of hydrodynamic simulations has been obtained. This is not meant to be a 
comprehensive discussion of final state interactions. It is only a description of one line of 
inquiry with which the author of this dissertation participated — namely the distorted wave 
emission function (DWEF) model. 



6.1 The RHIC HBT Puzzle 



Despite the success of early ideal hydrodynamical simulations of heavy ion collisions at 
RHIC, they had much difficulty fitting two-particle correlations while simultaneously match- 
ing single-particle data. For example, when the simulation parameters were set such that 
the experimental multiplicity and mean transverse momentum were at least roughly repro- 
duced, the predictions for Hanbury Brown and Twiss radii (recall [section 2.2.2P Rq and 
Rl are too large, while Rs is too small. In all reasonable cases, it seemed, the emission of 
particles occurred over a relatively long time period, causing the ratio Rq/Rs to be large. 
The experimental result, however, showed Rq/Rs ^ 1- 

This difficulty of describing HBT data with otherwise successful methods was dubbed the 
"RHIC HBT puzzle" |l4l|. It should be noted that — although adding viscosity improves the 
agreement of Rq /Rs in particular — even the most recent viscous hydrodynamic simulations 
do not completely resolve this puzzle j^] . Shedding light on this was the original motivation 
behind the development of the distorted wave emission function (DWEF) model, in which 



:ina^ 



142 



stat e interactions are introduced in the form of a relativistic optical rn odel formalism 



143]. Recent discussions of this puzzle can be found in Refs. |l44l . Il45 |. 
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6.2 The DWEF Model 

The medium produced in an ultrarelativistic heavy ion cohision is very hot and dense. It 
may be that the particles emitted from this opaque source experience significant refractive 
effects, which in turn affect measured quantities such as HBT radii in a way that cannot be 
captured by hydrodynamic models with only a simple Cooper-Frye freeze out mechanism. 
With this in mind, let us introduce a model for these final state interactions and develop 
the formalism for calculating various quantities within this model. Then by varying free 
parameters in the model one can show that it is possible to fit both single- and two-particle 
data for pions at RHIC, and then analyze the meaning of the values of the parameters 
necessary to fit the data. 

The DWEF formalism (along with many of the results) is laid out in detail in Ref. [142l |. 
Here we first present a general derivation of the formalism, in the manner of Ref. [J], and 
then we will describe the specific choices for the analytic form of the optical potential and 
emission function. 

6.2.1 Plane wave formalism 

For comparison, it is useful to start by deriving the formalism for calculating HBT radii in 
the absence of final state interactions, given some source function that represents particles 
that are emitted (incoherently) from the collision medium and which then travel without 
interaction into the detectors. In principle one could obtain this source function from some 
Cooper-Prye freeze out surface, but here it will be just a given function (and later an ana- 
ytic parametrization with tunable parameters). We follow one of the previous derivations 



146l | and then we can alter it appropriately to add an optical potential that the emitted par- 
ticles interact with. For simplicity we will specifically consider pions, the dominate hadron 
produced in a heavy ion collision. The extension to other particles is straight forward. 
We want to calculate the correlation function C(p, q) 

where P(pi, • • • Pn) is the probability of observing pions of momentum {pi} all in the same 
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event. The identical nature of ah pions of the same charge cause 



A state created by a random pion source \rj) is described by 



Ci 



P,P) = 2. 



147^ 



exp 



|0) = exp 



d-'pdt7]{p,t)-f{t)c^p)e''^^^ 



|0), (6.2) 



where t/j^ is the pion creation operator in the Heisenberg representation, 7(t) is the random 
phase factor that accounts for the chaotic nature of the source and c^(p) is the creation 
operator for a pion of momentum p. In particular, an average over collision events gives 



(7*(t)7(t')) 



6{t-t'), (6.3) 
6{h - h)6{t2 - U) + 6{h - U)5{t2 - H). (6.4) 



We note that as written, the state |r/) is not normalized to one. However, the normalization 
constant will divide out of the numerator and denominator of the correlation function. 
Therefore we do not make the normalization factor explicit here, but note that it enters 
when calculating the pion spectrum. 

For ip and its time derivative to obey the Heisenberg commutation relation one has 



^/e;e;\c{p),c\p')]=5^^\p-p'). 



Furthermore, we define 



d^xe-'P-''r]{x) 



(6.5) 



(6.6) 



The state \r]) is an eigenstate of the destruction operator in the Schroedinger representation, 
c(p): 



c{p)\v) 



En 



The correlation function is 



C(p,q) 



(7?|ct(p)ct(q)c(q)c(p)|T/) 



(6.7) 



(6.8) 



(r/|ct(p)c(p)|r?)(?7|ct(q)c(q)|r?) ' 
The use of Equation 6.4 and Equation 6.7| in the numerator of [Equation 6. 8| yields 

(77|ct(p)c^(q)c(q)c(p)|r?) = {v\c\p)c{p)\r]) {r]\c^q)c{q)\r]) + \{v\c\p)c{q)\rj)\'^ . (6.9) 
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Furthermore 



(77|ct(p)c(q)|ry) = J dt exp[-iiEp - Eq)t] 



V*{P,tHq,t) 

EpEq 

The quantity g{x,p) is denoted the emission function and is defined as 
g{x,p) = I dV7?*(x + ix',t)r/(x-ix',t)e^P-', 



so that 



(fn 1 1 

J^9{x, p)e-*P-^ = 7?*(x + -z, t)r/(x - -z, t) 

J (035((y + y')/2,t,p)e-^P-(y-^') =r,*(y,t)7?(y',t). 



(6.10) 

(6.11) 

(6.12) 
(6.13) 



The second expression appears in the right-hand-side of Equation 6.10| (if one uses Equation 6.6 ) 
so that we may write 



|c^(p)c(q)|??) 



exp[-i(p - q) ■ x] {p + q) ■ 



EpEq 



-9(x, 



(6.14) 



Using Equation 6.14 with p = q shows that the function g{x,p)/ Ep is the probabihty 
of emitting a pion of momentum p from a space-time point x. Using [Equation 6.9| and 



Equation 6.14 in Equation 6. 8| gives the desired expression: 



C(p q) - 1 / d'^x'g{x, ^K)g(x', ^K) exp[-ifc • {x - x') 

J d'^x d'^x'g{x,p)g{x' ,q) 



(6.15) 



where K = p + q and k = {Ep — Eq,p — q), and the factors of have canceled out. 

From a formal point of view, a key step in the algebra is the relation between the Heisen- 
berg representation pion creation operator i[}'^{x) and its momentum-space Schroedinger 
representation counterpart c^(p) that appears in Equation 6.2 



ip^{x) 

^(x) : 



d^p (p) 



d^p c(p) 



-ipx 



(27r)3/2 



(6.16) 
(6.17) 



(27r)3/2 

The operators c^(p) (c(p)) are coefficients of a plane wave expansion for ^/i^(x) (^/^(x)), with 
the plane wave functions e*P "^/(27r)3/2 being the complete set of basis functions. However, 
one could re-write ^''"(x) {i^{x)) as an expansion using any set of complete wave functions. 
We shall exploit this feature below. 
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6.2.2 Distorted waves — real potential 

We represent the random classical source, emitting pions that interact with a real, time- 
independent external potential U by the Lagrangian density: 

-C = 'ip^i-d'^ +U + m'^)'il^ + j{x)ip. (6.18) 



The current operator j{x) is closely related to the emission function g 147t |. In this La- 



grangian the terms U and jjx) are independent. Thus the relation between the emission 



function and U derived in 148l | need not be satisfied. Also note that the medium — and 
therefore also the potential — is in principle a time-dependent quantity. Nevertheless, for 
simplicity we take U to be time independent and it can be interpreted as a time-averaged 
quantity. 

The field operator tp^ can be expanded in the mode functions ipp ^ that satisfy: 

(-V2 +^/)4-)(x) =pVJ,"Hx). (6.19) 

These wave functions obey the usual completeness and orthogonality relations 

J d'p 4")* (x)4") (y ) = 5(3) (x - y) (6.20) 
J d3xV^{,-)*(x)V'^7^(x) = 6(^\p - p'), (6.21) 

so that one may use the field expansion 

i^{x)= y"d3p4-)(x,t)e-^^''*d(p), (6.22) 

with d^{p) being the creation operator for pions of momentum p in the basis of Equation 6.19 



These mode functions are the eponymous distorted waves which replace the plane waves 



of the previous section. The expansion Equation 6.22 assumes that U produces no bound 



states. If they did exist, the integral term would be augmented by a term involving a sum 
over discrete states. 

The availability of mode expansions when distortion effects are included means that the 
simplification of the correlation function can proceed as in the previous section. We again 
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use Equation 6.2 and Equation 6.4[ The use of the field expansion Equation 6.22 enables a 
generalization of the function ry(x): 



with 



so that 



\v) = exp 



d^xij^y ^(x,t)r/(x), 



|0). 



(6.23) 



(6.24) 



(6.25) 



Note that the ability to obtain a relation between the fj{p, t) and rj{x) rests on the relations 
Equation 6.20| and [Equation 6.21 



The state |ry) is an eigenstate of d{p). Thus the result 

Kr?|dt(p)d(q)|7?)p 



^^(p,q) = i + 



(r/|(it(p)d(p)|^)(ry|dt(q)d(q)|r?)' 



(6.26) 



very similar to Equation 6.8, is obtained. We need the matrix elements appearing in the 
numerator and find 



(r/|dt(p)d(q)|r?) 



^4^^3^/^^PL4i:^_M^(-)(x)V;(-)*(x')7?(x)r?(x',t). (6.27) 



EpEq 



and the use of Equation 6.13 allows us to obtain 



(r?|dt(p)d(q)|r?) 



1^ fdtd^xd^x'Ae^'(^^-^^^e-'P'-''' 



X 4 ^(x + x72)^^ ^*(x-x72)5(x,p'). 



(2^) 



(6.28) 



This result, which can be applied for p 7^ q and for p = q, specifies the evaluation of the 



correlation function of Equation 6.26 with the result 

\S{K,k)\' 



C(p,q) = l + 



S{p)Siq) 



(6.29) 



where 



S{K,k) = [ d^xdV-^e**(-^'-^p')e-^P'-^Vjr^(x + x72)V'i"^*(x-x72)5(x,p'),(6.30) 
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and 



S{p)^ I d4xdV|^e-^P'-''V(-)(x + x72)4--)*(x-x72)5(x,p') 



(6.31) 



This is the expression that is used in the DWEF formahsm jl42l . Il43l |. In principle one 
could use either Equation 6.15| or Equation 6.29 to analyze data, but the extracted space 
time properties of the source r/(x) would be different. 

A comment should be made on the possible momentum and energy dependence of the 
optical potential. The completeness and orthogonality relations are obtained with any 
Hermitian U which can therefore be momentum dep endent, but not energy dependent. 
As explained in section 5 (Equation 43) of Ref. [142], the real part of the potential can 
and should be thought of as a momentum-dependent, but energy-independent potential. If 
there were true energy dep endence a factor depending on the derivative of the potential 



with respect to energy 



149l | would enter into the orthogonality and completeness relations. 



6.2.3 Coupled channels 

If the optical potential U from the previous section is complex, the derivation above fails. 
Using the necessary completeness and orthogonality relations to relate r]{x) to f]{p, t) re- 
quires the use of a real potential. If we would like to include the effects of an imaginary 
part of the potential, we should investigate possible corrections to the above formalism. 

The optical potential or pion self-energy is an effective interaction between the pion 
and the medium. The medium is not an eigenstate of the Hamiltonian, but rather of Hq, 
which is the full Hamiltonian minus the Hermitian operator representing the pionic final 
state interactions. Eliminating the infinite number of possible states of Hq and representing 
these by a single state leads to a self-energy that is necessarily complex. Our procedure here 
is to specifically consider the infinite number of states of the medium, obtain a Lagrangian 
density that involves Hermitian interactions, and derive the optical potential formalism and 
any corrections to it. 

Let Pn denote a projection operator for the medium to be in a given eigenstate of Hq, 
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n. These obey 

n 

PnPm ~ ^n,mPn- (6.32) 

For the case of 7r-nuclear scattering, n would represent the nuclear eigenstates. Here n 
represents states of the medium in the absence of its interactions with pions. The correlation 
function is now given by 

where Pn(p) is the probability for emission of a pion of momentum p from the medium in 
a state n. Similarly Pn(p, q) is the probability for emission of a pair of pions of momentum 
p, q from the medium in a state n. The sums over n account for the inclusive nature of the 
process of interest. 

It is convenient to define the product of the field operator with the projection operator 

Pn- 

^n{x) = i^{x) Pn, (6.34) 

with 

V'(x) = ^V'n(x), (6.35) 

n 

using the complete nature of the set n. The Lagrangian density is given by 

- £ = ^ ^V^t . aV^n + 5^ V'i {{ml + M^)Snm + Unm) V-m + 3n{x)Mx), (6-36) 

n n,m n 

where 

^nm = ^mn = {^)nm (6.37) 

and U is the Hermitian interaction operator and M^, the m matrix element of the diagonal 
operator M^, represents the effects of the different energies of the states labeled by m. The 
field operator ip n can be expanded in the mode functions V'p,n • 

J2 W„^(x)47m(x, t) = (p2 + V2 - M2 - Unn{^))^i:l (x, t). (6.38) 
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Here the potential lA is taken as a local operator in the position space of the outgoing pion. 



To see the correspondence between the formulation of Equation 6.36 and Equation 6.38 



let V'l correspond to the field operator (and state) of the previous section and solve formally 
for ipp^m in terms of ifj^ I . It is convenient to define the operator lA with matrix elements 
given by 



Un,n' = (1 — ^n,l){'^ — ^n' ,l)Un,n' 



Then 



1 



(6.39) 



(6.40) 



where (V^ +p2 _ M'^)nm oc 6n,m, and is an operator giving when acting on the state 



n. Then rewrite Equation 6.38 in terms of ip^ ^ as 
The complex object 

1 



+ E ( 



V2 + p2 _ M2 - - ie 



/ nm 



V2 + p2 _ M 2 - - ie 



a non-local operator in coordinate space, can be identified with the optical potential, given 
by the operator V{Z) as a function of a complex variable Z: 

1 



V{Z) =Uii+ ^i™ 



v^ + z- ]vp -u 



(6.42) 



We proceed by employing Equation 6.35 and Equation 6.36 to compute the correlation 



function. The solutions of Equation 6.38 form a complete orthogonal set: 



E / (x)V'S(y) = 5(=^)(x-y) 

I d'x 4-^*(x)V'^7i(x) = 5(3)(p _ p'). 

The field expansion is now 

V'(X) = Jct'pY a(p)Pn47n(x) 



(6.43) 
(6.44) 

(6.45) 
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so that 



\v) = exp 



m) 



(6.46) 



where the state |0,m) is the pionic vacuum if the medium is in the state m, and rjn{x) 
represents the source for the state n. These state vectors obey the relations 



(0,n|0,m) = 5n,rn = (0,n|P„|0,m). 



(6.47) 



Define 



so that 



?/n(P,i)= / c^^a; r?„(x,f)^/;^,2*(x), 



exp 



d'pdt -,(t) t)aHp)P,.e"''' 



(6.48) 



(6.49) 



a(p)|r?)= [dtj{t)Y,^-^P^e^^r,t\^). 
The emission probabihty is given by 

n 

x472*(x)V'y(y) e^(^p-^«)* 



(6.50) 



(6.51) 



or using Equation 6.13 



EpE,{7^\a\p)a{cD\rj) = Y.n I dt d^x d^y J d'p' g„((x + y)/2, t, pOe^^P' ^-^) 

x472*(x)V'tn(y) e'iE^-^.)t, 



(6.52) 



where 



ff„(x,p) = y"(iV7?;(x+ ix',t)r/„(x- ix',t)e^P-^'. 



(6.53) 



If pionic final state interactions are ignored, the term enters and this may be identified 

with the emission function, g of previous sections. 
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The expression Equation 6.52 is the same as Equation 6.28 except that now we sum 
over the channels n. These sums may be expressed in terms of the optical model wave 
functions of Equation 6.40 , The term of [Equation 6.52 with n = 1 corresponds to the 



DWEF formalism, and the terms with n > 1 are corrections. As an example of a correction 
term, suppose part of the imaginary part of the optical potential arises from a pion-nucleon 
interaction that makes an intermediate A. Then a term corresponding to one of n > 1 
involves the emission of a pion from a nucleon that makes an intermediate A. 

It is difficult to assess the importance of the second term in a general way. The only 



obvious limit is that if states with n > 1 are not excited then Im{V) of Equation 6.42 must 
vanish. Conversely, if lm(y)=0, the states n > 1 must be above the threshold energy and 
the propagators that appear in the correction terms correspond to virtual propagation over 
a small distance with limited effect. 

Therefore in the following, results will be presented with the imaginary part of of the 
optical potential set at a vanishing value, and separately with any value allowed. While the 
former is correct within the model, the latter will have unknown corrections. Nevertheless, 
it will still be illustrative to look at both in the hope that it will give some idea of the effect 
of an imaginary part of the optical potential in addition to the more reliable information 
concerning the real part. 

6.2.4 Complete DWEF formalism 

The key unknown pieces in the expressions above are the emission function g and the optical 
potential Li. The wavefunctions ipp ^ can be calculated from Li, and then integrals can be 
performed to obtain the quantities of interest. Reiterating the results of lsection 6.2.2[ the 
correlation function for determining HBT radii is given by (recall [Equation 6.29 ) 



C(p,q) 



(/d4xd3x'|£;e-P'-'4-)(x + xV2)vi-^*(x-xV2)5o(x,pO) (p - q) ' 

(6.54) 

and single particle observables can be derived from one of the factors in the denominator: 
= y"d"xd3x'^e-^P'-"'4-'(x + x72)4-)*(x-x72)5o(x,p'). (6.55) 
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(In this sect ion, the notation So{x,p) = g{p,x) is used to make contact with the notation 



of Ref. 



We proceed by using an analytic parametrization t hat is inspired by hydrodynamic freeze 
out. A more detailed discussion can be found in Ref. 
The form used is: 



142] 



/ ^ o/ ^ cosh 7? -zpl 1 -(^~^o)' M± p{h) 

g(p,x) = So(x,p) = - — -ie2^ — e 2At^ —, ..J , (6.56) 



[7(b) = -{wQ + W2 p^) p{h). (6.57) 



p is the asymptotic pion momentum and M_|_ = yPj^ + rn^. Just as in previous chapters, 
it is natural to use Milne coordinates, although here we will use radial coordinates in the 
transverse plane: 



Tf] =arctanh(2;/t) r =\/ — b = vx^ + y^ 

(/) = arctan(y/x) h={b,(p). (6.58) 

Also as above, we will restrict ourselves to mid-rapidity data. 

p{h) represents the transverse density of the medium and is used for the transverse shape 
of both the emission function and the optical potential. It is normalized as p(0) = 1. The 
original DWEF model was restricted to rotationally symmetric systems (corresponding to 
central collisions) and used 

P(b) = p{b) = tViT^^ 7^- (6-59) 

This distribution has the correct exponential fall-off at large distance, and different choices 
of the parameters Rws and aws allow for a variety of shapes. To calculate the elliptic flow 
coefficient V2, it will be necessary to generalize this form for non-rotationally-symmetric 
systems (see [chapter 7 ). 



The velocity field that describes the dynamics of the expanding source in a central 
collision event is parametrized by a fluid rapidity r]t{h) 

u^{x) = {cosh T] cosh, rjt, sinhry^cos^, sinhr/t sin(/), sinhr/coshr/t). (6.60) 
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The flow rapidity is taken to have a Unear radial profile with strength rjf 

Vti^)=Vj^- (6-61) 

This also will have to be modified when calculating elliptic flow for a peripheral collision. 

The free parameters, then, are Arj, At, tq, fj.^, T, wq, W2, -RwSj ^ws, and r]f. These 
parameters were varied (with various of them occasionally held fixed) to rep rodu ce the 



single- and two-particle pion data for ^/s = 200 GeV Au-Au collisions at RHIC 142] 
6.2.5 Results for Central Collisions 



141 



143^ . With 



Calculations for central RHIC collisions were originally presented in Refs. 
the above insight, it is instructive to assess the possible importance of the imaginary part 
of the optical potential. 

A variety DWEF fits are performed, see ITable 6.T1 In Ref. [142l | the imaginary part of 



the optical potential as represented by the term W2 is about one tenth of the real potential. 
It is therefore possible that, in the limit that Im{w2) = 0, there would be no significant 
correction term, so we try to understand if removing the imaginary part of the optical 
potential can be done without degrading the quality of the fit. The results are shown in 



Figures [6.11 and [6121 An example of the previous calculations 142l . Il43l | is shown as the green 
dotted curve (second line of ITable 6.ip . The red solid curve (first line of lTable 6.ip shows the 
result of setting the imaginary potential to a vanishingly small value. This results in only 
a slightly worse description of the data. The changes in the imaginary part of the optical 
potential W2 are largely compensated by a reduction of the temperature from about 160 
MeV to about 120 MeV. We also point out that the length of the flux tube as represented 
by Ary is vastly increased, providing greater justification to our previous procedure of taking 
the length of the flux tube to be infinitely long in the longitudinal direction. However, the 
emissio n du ration is reduced to fm/c, which is similar to the results of the blast wave 



model |150| |. This means that all of the pionic emission occurs at a single proper time. 
This value justifies the use of a time-independent optical potential, but does seem to be 
difficult to understand because some spread of emission times is expected for a long-lived 
plasma. The results shown by the blue dashed curves (third line of lTable 6.1[) are obtained 
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Table 6.1: Four parameter sets obtained with slightly different procedures [4]. The values 
of represents the accuracy of the description of the data. 



T 




Ar 


Rws 


a-WS 




W2 


To 








(McV) 




(if) 


(fm) 


(fm) 


(fm-2) 








(MoV) 




121 


1.05 





11.7 


1.11 


0.495 


0.762 +0.0001i 


9.20 


70.7 


139.57 


300 


162 


1.22 


1.55 


11.9 


1.13 


0.488 


1.19+0.13i 


9.10 


1.68 


139.57 


117 


121 


1.04 


1.5 


11.7 


0.905 


0.564 


0.595 +0.0001i 


8.85 


70.7 


139.57 


451 


144 


0.990 


2.07 


12.57 


0.876 


0.0001 


0.0001+O.OOOH 


6.85 


oo 


83.5 


1068 



with fixing the emission duration to 1.5 fm/c, which is our previous value 142, Il43l |. The 
description of the spectrum is basically unchanged but the radii are less precisely described. 
The violet long-dashed curves (fourth line of lTable 6.ip show the DWEF fit using a vanishing 
optical potential. This does not give a good description of the momentum dependence of 
the radii and is associate with the largest deviation between our calculations and the data 
as represented by the values of lTable 6.11 

It is clear that the precision of our description of the data is improved by including 
the imaginary part of the optical potential. However, this is a quantitatively but not a 
qualitatively important effect. It is also true that including the real part of the optical 
potential is a qualitatively important effect. These results suggest that the correction terms 
embodied by the terms with n 7^ 1 of [Equation 6.52| are not very important, but non- 
negligible. It is also possible that an optical potential with a different geometry than the 
volume form that we have assumed might be able to account for the the neglected terms. 
However, an accurate assessment would require the development a theory that involves 
dealing with explicit models for gnijn and U. 
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Figure 6.1: Computed pionic spectrum. Red up right triangles: vr spectrum (STAR) Green 



inverted triangles: vr"'" spectrum points (STAR) 15ll ] Red solid line: DWEF fit with vanish- 
ing imaginary part of the optical potential, first line of lTable 6.1[ Green dotted line: DWEF 
fit including search on the imaginary part of the optical potential, second line of lTable 6.1[ 
Blue dashed line(almost entirely covered by the red solid curve): DWEF fit with vanishing 
imaginary part of the optical potential, Ar = 1.5 fm/c, third line of lTable 6.11 Violet long 
dashed line: DWEF fit including search on setting the optical potential to essentially 
0, fourth hue of lTable 6.ll 
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Figure 6.2: HBT radii. Curves are labeled as in 



Figure 6.1 



STAR data [152] 
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Chapter 7 
V2 IN THE DWEF MODEL 

Once results for central collisions have been calculated, the next interesting quantity is 
the elliptic flow coefficient V2 in non-central collisions. A few pieces need to be generalized 
for a non-rotationally symmetri c sys tem. For the transverse density p we take the modified 



Woods-Saxon profile from Ref. 



15C]. 



(exp[(-l)^] + l)2 
p(b) = ^'^^^ ^ , (7.1) 

(exp[(6J^ + ^-l)^] + l)^ 



with i?ws = yhi^x + ^y)- Thus lines of constant density in the transverse plane form 



ellipses with semimajor to semiminor axis ratio 

Next we must generalize the fluid velocity n, for which we again defer to Ref. | l5cj ]. 

n'^(x) = (cosh r/ cosh ?7f, sinh r]t cos (j)b, smh-qtsiiKph, sinhrycoshr/j). (7.2) 

The transverse direction is taken to be perpendicular to line s of constant density. It can be 



shown that the angle of such a fluid velocity, (j)^, obeys 



15C] 



0fe((A) =tan-i(--|tan,/.). (7.3) 

Ky 

The transverse fluid rapidity rit{h) is first taken to have the same elliptic symmetry 
as the density, increasing linearly with the "radial" coordinate b = — h I^ ™!*^) . 
Then added to this is a term proportional to cos(2(/)) representing the amount of elliptic 
fiow built up before freezeout 



' cos^ 6 sin^ 



r?i(b) =r/; 6J-^ + -^(l + a2 cos(2</))). (7.4) 



The momentum in these coordinates takes the form 



= (M_L cosh y,p_L cos sin M_L sinhy). (7.5) 
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Again we choose to focus on data at midrapidity, Y = 0, and so 



p ■ u = M± coshry coshryj — p± sinhr^j cos(0b — cpp). 



(7.6) 



Thus there are two extra parameters that characterize the departure from cyHndrical 
symmetry. In all, the parameters involved in this model are: Ar/, At, tq, /i^, T, wq, W2, Rx, 
Ryi flws, and 02. Rather than rerunning the fit for peripheral collisions (which would 
be prohibitively difficult numerically) we will here be interested in the effect of an optical 
potential like those found to give the best fit to central collision data, and therefore we will 
only adjust adjust and 02 to give reasonable results for non-central collisions. 

7.1 Calculating V2 

This section outlines how the calculation of V2 is carried out. A set of coupled differential 
equations must be solved numerically to obtain the wavefunctions tpp \ and then a five- 
dimensional integral must be performed (two of which can be done analytically with suitable 
approximations. ) 

7.1.1 The Wavefunctions ipp \x) 



ipp ^ satisfies Equation 6.19 Since U{h) is independent of t and z, we can write 



and Equation 6.19 becomes 



or 



Id 1 ^2 

+ 



bdb b'^dd)^ 



[/(b)+pi)v'i~)(b) = 0. 



Decomposing xpp ^ and f7(b) into angular moments 

00 

m=—oo 

C/(b)^Ef^"(^>"'^' 



(7.7) 

(7.8) 
(7.9) 

(7.10) 
(7.11) 
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results in 



db^^bdb h^^^^r"^ l^Unlm-nte 



^im4>^-im4,p = 0. (7.12) 



So the term in brackets vanishes identically for each m, and we must solve a set of 
coupled differential equations. In practice, every fm above a certain rrimax is set to zero, 
and a finite set of coupled equations is solved numerically. 

The boundary conditions are the same as for the cylindrically symmetric case — far out- 
side the medium one should have a canonically normalized plane wave plus an outgoing 
wave, i.e. 

fra{h » i?ws) = Jm{p b) + TmH^^ {p b) (7.13) 

with and Bessel functions and Hankel functions of the first kind, respectively. 

Details of this calculation can be found in appendix lA.il The program used to calculate 
the wavefunctions was tested in part by comparing to a semi-analytic solution described in 
appendix lA. 21 

7.1.2 Integration 

Once the wavefunctions are found, the integrals must be performed: 

„,.(c„s(2«) = i:*^f(^|fM. (7.14) 

jd(f)pS[p) 



with 



S{p) = I d^xd^x'^e-'^'-^'i^i'\^ + ^'/2)i;i->{^-^'/2)So{x,p') 



TdTd-qb db d(t) d^x'^^^e-'P'-'^'A'^ (x + x72)vi~^* (x - Ji.'/2)So{x, p') 



(7.15) 



Several approximations can make this more numerically tractable. If one assumes the 
the optical potential is approximately independent of the beam direction z as well as time, 
the T integral can be done analytically 



Tdre 2Ar^ = \/27rroAT. (7.16) 
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The p in tegral can also be done analytically with the following approximations (as in 
Ref. 142) 



1 cosh J] 



Jmax 



(7.17) 
(7.18) 



where the Bose-Einstein distribution is approximated by a sum over Boltzmann distributions 
truncated at some jmax , and so 

1 J 



/ 



dr] cosh rj e 



- cosh r]{ + cosh r\t ) 



2Ki 



At?- 



2 + ^M_L coshr/t 



(7.19) 



142l | 



Finally, we use the large source approximation 

V'i,;)(b + b72)47)*(b - b72) g(b'2) p. 47)(b)47)*(b) ^(b'^) exp (.Kx • b'), (7.20) 



with 



5(b'2) = 2 / d'K^ M_Lexp 



-M^ coshr/ 



exp [— zK^ • b'] . 



(7.21) 



After implementing all these approximations, for the numerator we have 



cos(2(?!)p) / d'^x S{p,x) 



2 ToM± 
{2-Kf 



E 



m,n,j 



X / dH pih)Uip,b) f*{p,b) e^^^'-^^^K^ 



1 , j 



2 + j;M± coshryt 



X J d(j)p cos{2(f)p)e' 



-i{m-n)(f>p ^j^pj^ sinh(r)t) cos((/)(,"(/)p) 



(7.22) 



and similarly for the denominator. The final three integrals are done numerically. 
More details of this part of the calculation can also be found in appendix lA.il 

7.2 Results 

We would like to determine the effect of adding final state interactions to hydrodynamic 
fits. To gain insight into this, we consider an emission function with parameter values taken 



from Refs 



.a 



142l |. which give the best description of the single particle data in general, 
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Table 7.1: Best fit parameter sets. The top line (Fit 1) is ageneral fit [142| | while the bottom 
line (Fit 2) is from a fit where Im{w2) is held at 0.0001 [4|. 





T 


Vf 


At 


Rws 




Wo 


W2 


-^0 


At] 




(MeV) 




i'-f) 


(fm) 


(fm) 


(fm-2) 




C-f) 




Fit 1: 


156.58 


1.310 


2.0731 


11.867 


1.277 


0.0693 


0.856+i0.116 


9.04 


1.047 


Fit 2: 


121 


1.05 





11.7 


1.11 


0.495 


0.762+iO.OOOl 


9.20 


70.7 



and also with the imaginary part of the optical potential held at zero (see Table 17.11 Also 
note that in both fits the chemical potential was fixed at the pion mass). 

We must make alterations to this central collision model to approximate a more pe- 
ripheral collision. The results for a central collision do not unambiguously imply what a 
peripheral collision will look like without appealing to a particular model for the dynamics 
of the system. We therefore choose reasonable parameters to approximately represent a 
collision with impact parameter ~ 7 fm, and then see how the resulting V2 depends on the 
strength of the optical potential. In principle one could vary all the parameters and do a 
separate fit of all the relevant experimental data (multiplicity, HBT radii, V2, etc.) for each 
of various collision centralities. However, the computing time to do so would be prohibitive, 
and here we are most interested in investigating only the effect of the interactions, so we 
proceed as follows. 

First, as in Ref. 



142] , we scale down i?ws > o-ws > and tq by the number of participants to 
the 1/3 power, with Npart taken from the Glauber model (with the same parameters used 
in Ref. [1]) for an impact parameter of and 7 fm {Npart = 377.5 and 171.544). Specifically 

Ft 

Rws 0.7688i?ws- Then we adjust the ratio such that the spatial eccentricity 



Rl + Rl 



(7.23) 



+ (a;2^ 

has a value of 0.035. This is a reasonable value corresponding to the spatial eccentric- 
ity at freezeout of hydrodynamic fits of peripheral collisions with impact parameter ~ 7 



fm. Note that the brackets in Equation 7.23 indicate a spatial average with weight given 



by Equation 7. 1 , while the spatial eccentricity in hydrodynamic simulations are typically 
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Figure 7.1: Calculated as a function of momentum with a2 = (a) and = 0.11,0.10 
(Fit 1, 2) (b). Points with err or ba rs are experimental data for pions at 20-30% centrality 



from the STAR Collaboration 



153|. 



given with respect to, e.g., energy density. We nevertheless keep the eccentricity from 



Equation 7.23 fixed at this value with an understanding that it is only a rough but still 



realistic guide to the shape. 

Lastly we must specify how much elliptic fluid flow is built up in earlier stages of the 
collision, represented by the value of 02 (recall [Equation 7.4] ) . First we set 02 = and see 
what V2 is generated by interactions with the optical potential in the absence of significant 
elliptic fluid flow ( [Figure T^TI ^a)). The calculated elliptic flow coefficient V2 is plotted as 
a function of momentum, along with the relevant experimental data. (Note that p in our 
calculation is the momentum of an asymptotically free pion detected far outside the medium, 
not the momentum of a particle as it is emitted inside the medium, and can therefore be 
compared directly to experiment.) Although we are only able to calculate up to a limited 
momentum, it is clear that flnal state interactions alone do not generate an appreciable 
value for V2 for either the general best-flt parameters (Fit 1) or those with a vanishing 
imaginary part of the optical potential (Fit 2). 

Next we increase 02 such that the experimental value for V2 is roughly obtained ( Figure 7T| b)). 
A value of 02 = 0.11 was required for the parameters from Fit 1, while 02 = 0.10 was suffi- 
cient to bring the emission function from Fit 2 into the physical regime. One can see that 
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the optical potential has a small but non-negligible effect — it decreases V2 on the order of 
10-25% of its zero-interaction value with a slightly smaller effect as momentum increases. 

7.3 Conclusion of DWEF V2 Calculation 

Final state interactions in the DWEF model were found to have a small, though not entirely 
insignificant effect on the elliptic flow coefficient V2- This is in addition to the indirect effect 
of adding flnal state interactions. For example, adding an optical potential changes other 
observables such as the multiplicity, which would alter parameters in a hydrodynamic flt 
such as freezeout temperature, which would then in turn have an effect on the calculated 
value of f 2 . 

The precise size of these effects in general can only be determined with a better under- 
standing of the model fits (e.g. Fit 1 versus Fit 2) in addition to a more detailed analysis — a 
full parameter search using all the relevant experimental data, or perhaps even by adding 
final state interactions directly i nto hydro dynamic simulations (i.e. a hydrodynamic after- 



burner in the vein of, e.g., Refs. 
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86, 



154( 1). It is reasonable, however, to conclude 
that final state interactions can affect the calculated value of V2 by as much as ~ 20% 
(in agreement with other investigations of final state interactions, e.g., Ref. |8J]), and so 
must be properly taken into account to have confidence in the quantitative predictions of 
hydrodynamic simulations at that level of precision. 



94 



Appendix 1 

DETAILS OF DWEF V2 CALCULATION 



A.l Numerical Implementation 



A program was written in C++, making use of the GNU Scientific Library (GSL) version 
1.9, to do the calculation of V2, as detailed here. 

The integral over the azimuthal angle of the pion momentum, (j)p is done as a sum 
using a simple trapezoid rule. This is because for each different value of (pp, a new set of 
differential equations must be solved. This also allows for the numerator and denominator 



of Equation 7.14 to be solved simultaneously, with just a factor of cos{2(pp) multiplied to 



the numerator when adding terms to the sum. 

For each term in the sum, then, first the wavefunctions Tpp ^ are obtained. They obey a 
set of coupled differential equations of the form 

S + II - ^ + ) ^™ - ^ ^" ^"^^'^'^ = ' ^^-'^ 

' n 

for all integers ra. This set is truncated, since large m moments (/m for m > p±Rws) con- 
tribute little to the wavefunction. Therefore, all fm for m greater than some rrimax are set to 
zero, leaving a finite {2mmax + 1) number of coupled ordinary differential equations. These 
are solved by calling a GSL solver. Using an embedded Runge-Kutta-Fehlberg method 



seemed to give the best performance. For these solutions. Equation 7.1 is integrated nu- 
merically to find the moments C/„. This is done with the GSL adaptive integration routine 
for oscillatory functions. 

To match to the proper boundary conditions, one must find {2mmax + 1) linearly inde- 
pendent solutions to this set of equations and take the correct linear combination of these 
solutions that matches the desired boundary conditions. The straightforward choice for 
these linearly independent solutions is to sequentially solve for the case where only one of 
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the partial waves is non-zero near the origin. For example, for the n'th solution let: 

P 

and then solve the set of differential equations up to some cirbitrcirily large hfnax 

far outside 

the potential. We can then match each partial wave in this n*^ solution to the form: 

) = Am,nJmip b) + Bm,nHg\p b). (A.3) 

The final wavefunction is then given by the linear combination of these solutions that 



matches the form of Equation 7.13 at bmax'- 

fm{b) = Y.CnfmAb)- (A.4) 

n 

This part of the program was tested with the trivial case of zero optical potential, in 
addition to comparing to a separately written program that calculates only the cylindrically 
symmetric well as to the results of the semi-analytical test case described in appendix 

Once these wavefunctions are obtained and stored in memory, the integral over h and 
(p in Equation 7.22 can be performed in addition to the sum over Boltzmann factors. The 
integrations are done with two GSL adaptive integration routines, one embedded in the 
other. The sum is done inside the argument of the integrals. 

A. 2 Semi- Analytic Test Case 

To test the numerics, the case of a pion moving through an elliptically-shaped step-function 
potential was solved (semi-) analytically making use of elliptic coordinates. This can be 
compared to the case of a^s (see [chapter 7 ) . 



We want to solve [Equation 6.19| with U{h) an elliptically shaped step function — a finite 
potential inside an ellipse in the transverse plane, with zero potential outside. 

It is useful to change to elliptic (cylindrical) coordinates, denoted u and v. Think of u 
as a 'radial' coordinate that runs from to oo and v as an 'angular' coordinate that runs 
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from to 2tt 



X = a cosh(ti) cos{v) 
y = a sinh(M) sin(t>). 



(A.5) 



Note the major and minor axes of the resulting confocal elhpses are reversed from the shape 
of the density used in the main calculation (which is larger in the y direction). This is to 
maintain consistency with the conventional definition of elliptic coordinates. At the end 
one can simply take (pp {(pp + vr) to match the usual convention in RHIC papers. 
Consider the case 

U{h) = U{u) = Uo e(no - u). (A.6) 
The sharp boundary at u = uq is an ellipse with major and minor axes 

Rx = a cosh(tio) 

Ry = a sinh(uo). (A- 7) 

In this coordinate system the Laplacian is 



V 



1 



(sinh^(ti) + sin^(u)) \du'^ dv"^ J 



and so [Equation 6.19 becomes 

1 



(sinh^(u) + sin^(T;)) \du? dv'^ J 



U{u) +p2 



V'p(b) = 



or equivalently 



with 



^2 

-^ + 2q{u) cosh(2u) + — - 2g('u) cos(2t>) 



V'p(b) = 



An) 



- U{u)) . 



(A.9) 



(A.IO) 



(A.ll) 



On the inside of the potential and on the outside separately, q{u) does not depend on u and 
these cases can be solved with separation of variables and the solutions patched together at 
u = Uq. Let 



2 

Qout = -jP ■ 



(A.12) 
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Start by expanding ipp{h) in terms of so-called ellip tic s ines and cosines of the 'angular 



variable v. They are solutions of 'Mathieu's equation' 



155|: 



„ + 2q cos(27;) C(a, q,v) = a C(a, q, v). (A.13) 
au^ / 

The general solutions are called 'Mathieu functions,' usually denoted C(q, q, v) for solutions 
even in the coordinate v and S{a^ q, v) for odd. Demanding periodicity of the variable v 
allows only certain discreet eigenvalues a (denoted here an for the even functions and /?„ 
for the odd functions). This (complete) set of periodic solutions is commonly called elliptic 
sines and elliptic cosines: 



C{an,q,v) = cen{v,q) 
S{(3n,q,v) = sen{v,q). 



(A.14) 



The general solution of Equation A. 10 can be written in terms of these elliptic sines and 
cosines: 



V'p(b) = ^ [fcMcen{v,q) + fsMsen{v,q)] 



n=0 



Plugging this in to Equation A. 10 gives 



dv? 



+ 2q cosh(2n) — ar, 



+ 2q cosh(2ti) — 



/c„(n)=0 

/.J^) = o. 



(A.15) 

(A.16) 
(A.17) 



This is called the modified Mathieu equation, which can be obtained from Equation A.13| 
by replacing v ^ {i u). Note that the eigenvalues are different for the functions correspond- 
ing to ce„ and se„ (/c„ and /^^ above, respectively). The general solution is then the same as 
for the original Mathieu equation, analytically continued with v ^ (i u), though typically 
they are organized by boundary co nditi ons analogous to Bessel and Neumann functions 
(denoted Jen{u,q), Nen{u,q), etc.) [156l |: 



fcM = Cc„Jen{u,q) + Sc„Nen{u,q) 
fsM = Cs„Jon{u,q) + Ss„Non{u,q). 



(A.18) 
(A.19) 
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Note that there are many different sets of so-called Mathieu functions, each being a 
complete orthogonal basis. Replacing gj„ with qout results in a different basis, and there are 
separate sets of modified Mathieu functions corresponding to the eigenvalues of the elliptic 
sines and elliptic cosines (a„ and 

By requiring continuity at the u = line segment one finds that the general solution 
inside the potential is: 

v)=Y^ [Ce^Jeniu, qin)cen{v, qin) + Co^Joniu, qin)sen{v, qin)] (A.20) 

n 

with undetermined coefficients Ce^",Co^". 



157| 



Outside, we write the solution as the sum of a plane wave and an outgoing wave 

^[( — Jen{u,qout) + Ce'^^He^^\u,qout)] cen{v, qout) cen{4>p, Qout) 

+ (^-^Jon{u,qout) + CoT'Ho^^^ iu,qout)^ sen{v,qout)sen{4>p,qout)], 

where the H's are analogous to Hankel functions 

He^^\u,q) = Jen{u,q)+i Nen{u,q) (A.22) 
Ho^^^ {u, q) = Jon{u, q) + i Non{u, q) (A.23) 

and the plane wave coefficients pn and s„ are 

dv e'P''cen{v,qout) (A.24) 

/ dv e'P''seniv,qout)- (A.25) 
Jo 

The coefficients Ce°'"* and Cd^^, along with the analogous 'inside' coefficients are deter- 
mined by matching boundary conditions. 

To match at the u = uq boundary, project the 'inside' angular functions 
(e.g. cen{v,qin)) in terms of the 'outside' ones (e.g. cen{v , qout)) ■ 

oo 

cej{v, qin) = ^ Bj^cen{v, qout) (A.26) 



1 


1 f^TT 


Pn 


Jo 


1 






71" Jo 



n=0 
oo 



sej{v, qin) = ^ Bj^sen{v, qout), (A.27) 



n=0 



99 



with 

1 Z"^'' 

^jn = - dvcej{v,qin) cen{v,qout) (A.28) 

^ JO 

1 r^'' 

^jn=- dvsej{v,qin) sen{v,qout)- (A.29) 
Then the 'inside' wave functions are 



^p" = XI i^^T -^^J^"' ^«") ^jncen{v, qout) + Cof Joj{u, qin) B]„sen{v, qout)] ■ (A.30) 

The coefficients (Ce™, Co™, Ce°"*, Cof^^) can then be determined by demanding that ■i/' 
and its gradient be continuous at u = uq, which gives the fohowing relations: 



^ CefJej{uQ, qin)Bj„ = —Jen{uo, qout)cen{(j)p, qout) + Ce^^He^^^ {uq, qout)cen{(j)p, qout) 

(A.31) 

^ CofJojiuQ, qin)Bjn = —JoniuQ, qout)sen{<t>p, qout) + Co°^*Ho''^\uo, qout)sen{(t>p, qout) 

(A.32) 

^ CefJe'jiuQ, qin)Bjn = —Je'^{uo, qout)cen{^p, qout) + Ce'^*He'^^\uo, qout)cen{<t>p, qout) 

(A.33) 

'^CoJ'JOj{uo,qin)Bjn = —Jo'^{uo, qout) sen{(l)p, qout) + Co^*Ho'^^\uo,qout)sen{(l)p,qout)- 

(A.34) 

The plane wave coefficients {pn, Sn) as well as the coefficients from the projection (-Bj„, 
Bj^) must be solved numerically. In addition, to compare to the fm in the main calculation, 
the resulting wavefunctions are integrated to project out the usual angular moments. Hence 
the description as a "semi-analytical" test case. In fact, this implementation (done in 
Mathematica) saves no time over the original numerical version, but it does provide an 
independent check. 
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Appendix 2 

NOTATION, CONVENTIONS, AND DEFINITIONS 

All notational definitions axe defined when first introduced, but frequently used notation 
is collected here for easy reference (or at least notation that is used in well-separated parts 
of the manuscript). 

• All quantities are reported using a system of units such that c = h = ks = 1 { "natural 
units"). I.e., all velocities are measured as fractions of the speed of light c, etc. 

• The space-time metric in flat space is taken as 5^,^ = diag(l, -1, -1, -1), such that 
timelike 4-vectors have positive norm and spacelike vectors negative. 

• Projectors: 

AM'^ = ^'^'^ - (B.l) 

p^; ^ a^a;^ + A^Aj; - ^a^-a,^ , (b.2) 

such that Ufj^A^^'^ = u^P^'^ = gixvP^p = 0- Projecting with A'*'' makes a quantity 
transverse to a fluid velocity u^, while P^^ makes it transverse, traceless, and sym- 
metric under interchange of indices. 

• Derivatives: 

D = u^'df, , (B.3) 
V^^A^Sa, (B.4) 

so that dn = u^D + V^^. In the fluid rest frame, these are the time derivative and 
spatial gradient, respectively. I.e., in the non-relativistic limit 

DKdt + v-d + 0{\v\'^) , (B.5) 
V^-d + 0{\v\). (B.6) 
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• Brackets: 

^ 1 + A^B""^ , (B.7) 

^["5/3] ^ 1 - , (B.8) 

^(a^/3> ^ p^Pa^'B'' , (B.9) 

which are used to define cr'^'^ = V^^^u"^ and the fiuid vorticity oj'^^ = — Vl'^u^'l. 
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